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ABSTRACT 

Binary neutron star (BNS) mergers are the leading model to explain the phenomenology of short gamma-ray 
bursts (SGRBs), which are among the most luminous explosions in the universe. Recent observations of long- 
lasting X-ray afterglows of SGRBs challenge standard paradigms and indicate that in a large fraction of events a 
long-lived neutron star (NS) may be formed rather than a black hole. Understanding the mechanisms underlying 
these afterglows is necessary in order to address the open questions concerning the nature of SGRB central 
engines. However, recent theoretical progress has been hampered by the fact that the timescales of interest 
for the afterglow emission are inaccessible to numerical relativity simulations. Here we present a detailed 
model to bridge the gap between numerical simulations of the merger process and the relevant timescales 
for the afterglows, assuming that the merger results in a long-lived NS. This model is formulated in terms 
of a set of coupled differential equations that follow the evolution of the post-merger system and predict its 
electromagnetic (EM) emission in a self-consistent way, starting from initial data that can be extracted from 
BNS merger simulations and taking into account the most relevant radiative processes. Moreover, the model 
can accomodate the collapse of the remnant NS at any time during the evolution as well as different scenarios 
for the prompt SGRB emission. A second major reason of interest for BNS mergers is that they are considered 
the most promising source of gravitational waves (GWs) for detection with the advanced ground-based detector 
network LIGO/Virgo coming online this year. Multimessenger astronomy with joint EM and GW observations 
of the merger and post-merger phase can greatly enhance the scientific output of either type of observation. 
However, the actual benefit depends on whether a suitable EM counterpart signal to the GW emission can be 
identified (ideally bright, isotropic, long-lasting, and associated with a high fraction of BNS merger events). 
The model presented here allows us to search for such counterparts, which carry the signature of a long-lived 
remnant NS. The present paper is devoted to a detailed discussion of the formulation and implementation of the 
model. In a companion paper, we employ this model to predict the EM emission from ^ 10“^ to ^ 10^ s after 
a BNS merger, considering a wide range of physical parameters, and discuss the implications in the context of 
SGRBs and multimessenger astronomy. 

Keywords: gamma-ray burst: general — gravitational waves — pulsars: general — radiation mechanisms: 
general — stars: magnetars — stars: neutron 


1. INTRODUCTION 

The coalescence of binary neutron stars (BNS) represents 
the most promising source of gravitational waves (GWs) 
for the detection with ground-based interferometric detectors 
such as advanced LIGO and Virgo (Harry et al. 2010; Acca- 
dia et al. 2011). At the same time, BNS mergers are respon¬ 
sible for observable electromagnetic (EM) emission that car¬ 
ries complementary information on the source and that can 
be combined with the GW detection to significantly enlarge 
the scientific output. As an independent observational chan¬ 
nel, EM signals provide positional and temporal information 
that can enhance the search sensitivity of the GW detectors 
and lead to GW detections (e.g., Abadie et al. 2012b; Aasi 
et al. 2014; Williamson et al. 2014; Clark et al. 2014). In¬ 
versely, EM follow-up observations triggered by a GW de¬ 
tection can confirm the astrophysical origin of the event (e.g., 
Evans et al. 2012; Abadie et al. 2012a; Singer et al. 2014), 
provided a suitable, characteristic EM counterpart signal can 
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be identified with the merger and/or post-merger phase of the 
evolution of the BNS system. In addition, EM counterparts 
can constrain physical properties of the BNS merger remnant 
and its dynamics that cannot be probed by GW observations 
(e.g., magnetic fields, mass ejection). In the case of a confi¬ 
dent GW detection, important contraints on the source prop¬ 
erties, the EM emission mechanisms, and the energetics in¬ 
volved can be obtained even from a non-detection of an EM 
counterpart. Moreover, detecting such EM signals could sig¬ 
nificantly improve the sky localization of the source and thus 
the identification of the host galaxy, allowing for two indepen¬ 
dent measurements of the redshift and the Hubble constant 
(e.g., Schutz 1986; Metzger & Berger 2012; Berger 2014). 
Eurthermore, joint EM and GW observations can reveal im¬ 
portant information on when and how short gamma-ray bursts 
(SGRBs) can be produced in BNS mergers. In particular, they 
can verify or falsify the recently proposed ‘time-reversal’ sce¬ 
nario for SGRBs and possibly provide a very accurate method 
to determine the lifetime of a remnant NS and to place strong 
constraints on the unknown equation of state of nuclear matter 
at high densities (Ciolfi & Siegel 2015b, a). With the advanced 
LIGO/Virgo detector network starting their first science runs 
later this year, such multimessenger astronomy will become 
reality in the very near future. The actual benefit of joint GW 
and EM observations, however, depends on our knowledge of 
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the expected EM signals associated with BNS mergers and 
the underlying physical mechanisms, which are still poorly 
understood. Therefore, a deeper understanding of these EM 
signals is urgently needed. 

BNS mergers are commonly considered the leading sce¬ 
nario to explain the phenomenology of SGRBs (e.g., Paczyn- 
ski 1986; Eichler et al. 1989; Narayan et al. 1992; Barthelmy 
et al. 2005; Eox et al. 2005; Gehrels et al. 2005; Shibata et al. 
2006; Rezzolla et al. 2011; Paschalidis et al. 2015). Currently, 
very strong evidence for the association of SGRBs with BNS 
or neutron star-black hole (NS-BH) mergers comes from the 
recent detection of possible radioactively powered kilonova 
events (e.g., Li & Paczyhski 1998; Kulkami 2005; Rosswog 
2005; Metzger et al. 2010; Tanvir et al. 2013; Berger et al. 
2013; Yang et al. 2015). While this is still a matter of de¬ 
bate, joint EM and GW observations will provide a power¬ 
ful tool to unambiguously confirm the association of SGRBs 
with compact object coalescence. Eor this reason, the prompt 
7 -ray emission of SGRBs lasting less than « 2 s has so far 
been the prime target of studies on joint EM and GW obser¬ 
vations (e.g., Evans et al. 2012; Abadie et al. 2012a; Singer 
et al. 2014; Abadie et al. 2012b; Aasi et al. 2014; Williamson 
et al. 2014; Clark et al. 2014). 

Nevertheless, the SGRB prompt emission is thought to be 
collimated and will thus be beamed away from the observer 
in most cases (see Berger 2014 for an overview of observa¬ 
tions to date; see also Section 7 of Siegel & Ciolfi 2015a). 
During many years of observations, including the 10 years 
of operation of the Swift satellite (Gehrels et al. 2004), no 
SGRB with known redshift has been detected within the sen¬ 
sitivity volume of advanced LlGOWirgo. Eurthermore, joint 
EM and GW observations are so far based on the assumption 
that the EM and the GW signals are characterized by a rela¬ 
tive time lag of at most a few seconds. However, the details of 
how the prompt emission in SGRBs is generated still remain 
unclear, and the burst could, e.g., be generated a long time 
(~ 10^ — 10^ s) after the merger (Ciolfi & Siegel 2015b, a). 
In the latter case, joint observations focusing on a time win¬ 
dow of a few seconds around the time of merger would miss 
the SGRB. Likewise, focusing on a time lag of a few sec¬ 
onds around a detected SGRB could lead to a non-detection 
of GW emission. Hence, when assuming coincidence within 
a short time window around the time of merger, one has to 
be cautious when drawing astrophysical conclusions in either 
case. Einally, even if a SGRB is observed in coincidence with 
a GW signal, such an observation alone will unlikely be able 
to distinguish between a BNS and a NS-BH merger, as the 
SGRB emission is expected to be very similar for both pro¬ 
genitor models (see below). It is therefore important to iden¬ 
tify bright, long-lasting and highly isotropic EM counterparts 
that are produced in a high fraction of events and that can dis¬ 
tinguish between a BNS and a NS-BH progenitor system. 

A standard model to explain the generation of the SGRB 
prompt emission in BNS and NS-BH mergers is an accre¬ 
tion powered relativistic jet from a BH-toms system that is 
formed soon (< 10 — 100 ms) after merger (e.g., Narayan 
et al. 1992; Shibata et al. 2006; Rezzolla et al. 2011; Pascha¬ 
lidis et al. 2015). This accretion process and the resulting 
energy release cease once the torus has been accreted on a 
timescale of at most one second, which is consistent with 
the typical timescale of the prompt SGRB emission (< 2 s). 
However, recent observations by the Swift satellite have re¬ 
vealed long-lasting X-ray afterglows in a large fraction of 
SGRB events that are indicative of ongoing energy ejection 


on much longer timescales up to ~ 10^ s (Rowlinson et al. 
2013; Gompertz et al. 2013, 2014; Lii et al. 2015). Even 
considering that the interaction of the jet with the interstellar 
medium might produce an afterglow signal lasting longer than 
the accretion timescale, persistent emission with a duration of 
^ 10^ s remains difficult to explain within this BH-torus sce¬ 
nario (Kumar & Zhang 2015 and referecnes therein). A pos¬ 
sible alternative is that a large fraction of BNS mergers lead 
to the formation of a stable or at least sufficiently long-lived 
NS rather than a BH-torus system (e.g., Zhang & Meszaros 
2001; Metzger et al. 2008; Rowlinson et al. 2010; Bucciantini 
et al. 2012; Rowlinson et al. 2013; Gompertz et al. 2013; Lii 
et al. 2015). Such a long-lived NS can power ongoing en¬ 
ergy ejection on the relevant timescales via loss of rotational 
energy. This is a clearly distinctive feature of BNS mergers 
as opposed to NS-BH mergers, which cannot produce a rem¬ 
nant NS. NS-BH mergers are thus challenged as a progenitor 
model for at least a large class of SGRBs. 

The formation of a long-lived NS is indeed a very likely 
outcome of a BNS merger. The merger product depends on 
the masses of the two progenitor NSs and the equation of state 
of nuclear matter at high densities, which is unknown. Recent 
observations of high-mass NSs (Demorest et al. 2010; Anto- 
niadis et al. 2013) indicate a maximum gravitational mass for 
stable NSs of Mtov 2M0. However, NSs with masses 
M > Mjoy can be centrifugally supported against gravita¬ 
tional collapse by uniform rotation up to the mass-shedding 
limit, which is known as the supramassive regime. This re¬ 
sults in a maximum mass for uniformly rotating configura¬ 
tions of Mjupra ~ 1.2Mxov ^ 2.4M0 (Lasota et al. 1996). 
Eurthermore, the distribution of NSs in binary systems is 
sharply peaked around 1.3 — 1.4M0, with the first bom NS 
slightly more massive than the second one; this leads to a typ¬ 
ical remnant NS mass of « 2.3 — 2.4 M 0 when accounting for 
neutrino losses and mass ejection (Belczynski et al. 2008). 
Progenitors of lighter NSs are much more abundant than pro¬ 
genitors for massive NSs, and population synthesis calcula¬ 
tions show that 99% of all BNS mergers should lead to a rem¬ 
nant mass between « 2.2 — 2.5 M 0 (Belczynski et al. 2008). 
Hence, the most likely product of a BNS merger should be 
a supramassive NS. Depending on Mtov there might also be 
a significant fraction of stable NSs. Eurthermore, some BNS 
mergers would lead to a slightly hypermassive NS (i.e., above 
the maximum mass supported by uniform rotation), which, 
however, could still migrate to a supramassive configuration 
through subsequent mass loss (cf. Section 4.1). Only a small 
fraction of BNS mergers should promptly form a BH-torus 
system as assumed in the standard model. 

EM emission from a long-lived remnant NS when applied 
to a large class of X-ray afterglows of SGRBs has so far been 
modeled in a simple way as dipole spin-down emission from 
a uniformly rotating magnetar (e.g., Rowlinson et al. 2013; 
Gompertz et al. 2013; Lii et al. 2015). This model assumes 
an instantaneous and direct conversion of spin-down lumi¬ 
nosity Lsd into observed X-ray luminosity Lx by some un¬ 
specified process, oc. Lx, and consists of a simple an¬ 
alytically specified formula to fit the X-ray lightcurves. In 
particular, it does not take into account baryon pollution due 
to dynamical mass ejection and subsequent neutrino and mag¬ 
netically driven winds (e.g., Hotokezaka et al. 2013; Oechslin 
et al. 2007; Bauswein et al. 2013; Kastaun & Galeazzi 2015; 
Dessart et al. 2009; Siegel et al. 2014; Metzger & Eemandez 
2014), which leads to a much more complex post-merger evo¬ 
lution (see, however, Metzger et al. 2008; Bucciantini et al. 


Electromagnetic emission erom long-lived BNS merger remnants I 


3 


2012). Based on a simple dynamical model, Yu et al. (2013) 
and Gao et al. (2015) have investigated EM emission from a 
long-lived millisecond magnetar surrounded by an envelope 
of previously ejected matter, finding a late-time brightening 
(termed “magnetar-driven macro-nova” ) consistent with the 
optical and X-ray (re)brightening of GRB 080503. Based on 
a refined physical model, Metzger & Piro (2014) have investi¬ 
gated the evolution of a stable millisecond magnetar in a sim¬ 
ilar setup and reported late-time brightenings of the luminos¬ 
ity in the optical and X-ray band, compatible with features 
observed in GRB 080503 and GRB 130603B (see Section 6 
for a comparison of these previous models to ours). Accu¬ 
rately modeling the EM emission from a long-lived remnant 
NS is challenging because of the complex physics involved 
(baryon pollution, thermal effects, neutrino emission, strong 
magnetic fields with complex structure, strong differential ro¬ 
tation) and the wide range of timescales involved (from ms 
to days after the merger). While the early post-merger phase 
and thus the generation of the prompt SGRB emission can 
be probed by numerical relativity simulations, such typical 
timescales for EM afterglows are inaccessible to those simu¬ 
lations. On the other hand, semi-analytical modeling has so 
far concentrated on computing EM emission in the late-time 
regime (~ 10^ — 10® s after merger; Yu et al. 2013; Gao et al. 
2013; Metzger & Piro 2014; Gao et al. 2015). 

Here and in a companion Paper (Siegel & Ciolfi 2015a, 
henceforth Paper II), we consider the likely case of a long- 
lived NS remnant and provide a dynamical model to self- 
consistently compute the EM emission of the post-merger sys¬ 
tem based on some initial data that can be extracted from 
a numerical relativity simulation tens of milliseconds after 
the merger. This model thus bridges the gap between the 
short timescales accessible to numerical relativity simulations 
and the timescales of interest for EM afterglow radiation as 
recorded by satellite missions like Swift. We note that our 
model does not include the EM emission associated with the 
formation of a relativistic jet, which can be added to the EM 
signal predicted here. The model is very general and should 
be applicable to any BNS merger that leads to the formation of 
a long-lived NS. As we have argued above, this should cover 
the vast majority of BNS merger events. Therefore, our model 
represents an important tool to study and identify promising 
EM counterparts of BNS mergers for coincident EM and GW 
observations as discussed above and to investigate the nature 
of long-lasting X-ray afterglows observed in a large fraction 
of SGRB events. The present paper is devoted to a detailed 
discussion of the physical model and its numerical implemen¬ 
tation. In Paper II, we apply the model to a large number of 
possible long-lived BNS merger remnants, employing a wide 
range of physical input parameters. Our results and their as- 
trophysical implications are discussed in Paper II. 

This paper is organized as follows. Section 2 describes 
the phenomenology underlying our evolution model. In Sec¬ 
tion 3, we provide a brief summary of the model, which is 
formulated in terms of sets of highly coupled ordinary differ¬ 
ential equations. The subsequent section describes the ingre¬ 
dients to these equations in detail. In Section 5, we discuss 
some numerical aspects for integrating the model equations 
and, in particular, present a scheme to reconstruct the observer 
lightcurves and spectra including relativistic beaming, the rel¬ 
ativistic Doppler effect and the time-of-flight effect. Section 6 
is devoted to discussion and conclusions. Several appendices 
are added to streamline the discussion in the main part of the 
paper. 


2. PHENOMENOLOGY 

This section is intended to provide a conceptual outline of 
the phenomenology our evolution model is built upon. We as¬ 
sume that a BNS merger leads to the formation of a long-lived 
NS and divide the post-merger evolution into three phases: an 
early baryonic-wind phase (Phase I; f 1 — 10 s), a pulsar 
wind shock phase (Phase 11; f ^ 10 s), and a pulsar wind 
nebula (PWN) phase (Phase III; ~ 10 s < f < 10^ s). These 
evolution phases are depicted in Eigure 1 . Our model provides 
a self-consistent dynamical evolution once the physical prop¬ 
erties of the system in the early post-merger phase are speci¬ 
fied. These properties can be extracted or estimated from BNS 
merger simulations in general relativity. In particular, we start 
the evolution at a few to tens of milliseconds after the BNS 
merger, once a roughly axisymmetric state of the remnant NS 
has been reached and the strong GW emission characterizing 
the merger phase has been severely damped. 

In our reference scenario. Phase I starts with a supramassive 
NS, which is a very likely outcome of a BNS merger (see Sec¬ 
tion 1). Nevertheless, the phenomenology described here also 
applies if the NS is stable or if its mass is only slightly above 
the hypermassive limit, provided that in the latter case the star 
enters the supramassive regime through additional mass loss 
before collapsing to a black hole (see below). 

The newly-bom NS is differentially rotating and endowed 
with strong magnetic fields (up to magnetar field strengths) 
due to compression of the two progenitor stars and magnetic 
field amplification mechanisms acting during the merger, such 
as the Kelvin-Helmholtz instability (e.g.. Price & Rosswog 
2006; Zrake & MacFadyen 2013; Giacomazzo et al. 2015). 
Moreover, further amplification occurs during Phase I, i.e., 
as long as differential rotation is active, via magnetic wind¬ 
ing and possibly the magnetorotational instability (e.g., Duez 
et al. 2006; Siegel et al. 2013; Kiuchi et al. 2014). These latter 
mechanisms are likely to dominate the removal of differential 
rotation itself, which occurs on the timescale fdr- In this case, 
fdr is roughly given by the Alfen timescale and can be as long 
as fdi 1 — 10 s for the objects considered here, depending 
on the initial magnetic field strength (Shapiro 2000). 

At birth, the remnant star is already surrounded by material 
dynamically ejected during or shortly after the merger, with a 
total mass of up to Mgj < 10“® Mq (e.g., Hotokezaka et al. 
2013; Bauswein et al. 2013; Kastaun & Galeazzi 2015; see the 
discussion in Section 4.1). However, additional mass ejection 
is expected to take place during Phase I, which can even domi¬ 
nate over this early dynamical outflow in terms of total ejected 
mass. One main mass ejection mechanism results from mag¬ 
netic field amplification in the stellar interior, which causes a 
build-up of magnetic pressure in the outer layers of the star. 
This pressure rapidly overcomes the gravitational binding at 
the stellar surface, launching a strong baryon-loaded magne¬ 
tized wind (Siegel et al. 2014; Siegel & Ciolfi 2015c). Fur¬ 
thermore, substantial mass loss can be caused by neutrino- 
induced winds over typical timescales for neutrino cooling of 
^ 1 s (Dessart et al. 2009). For both mechanisms, i.e., 
magnetically and neutrino-induced winds, we expect highly 
isotropic mass ejection. The material ejected from the NS sur¬ 
face is typically hot with temperatures of up to tens of MeV 
(Siegel et al. 2014; Kastaun & Galeazzi 2015), while further 
energy is carried by the wind in the form of a strong Poynt- 
ing flux (Siegel et al. 2014). These winds transport material 
outward at speeds of at most Uej.in ^ 0.1c, creating a hot and 
optically thick environment. EM emission from these radi- 
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Figure 1. Evolution of the system according to the proposed scenario (with 
increasing spatial scale). A BNS merger (top left) forms a differentially ro¬ 
tating NS that emits a baryon-loaded wind (Phase I). The NS eventually set¬ 
tles down to uniform rotation and inflates a pulsar wind nebula (or simply 
‘nebula’) that sweeps up all the ejecta material into a thin shell (Phase II). 
Spin-down emission from the NS continues while the nebula and the ejecta 
shell keep expanding (Phase III). 


ally expanding winds is expected to be predominantly ther¬ 
mal, due to the very high optical depths at these early times. 
However, because of the high optical depth, radiative energy 
loss is still rather inefficient. 

As differential rotation is being removed on the timescale 
fdr. the NS settles down to uniform rotation. Mass loss is 
suppressed and while the ejected matter keeps moving out¬ 
ward the density in the vicinity of the NS is expected to 
drop on roughly the same timescale. In the resulting essen¬ 
tially baryon-free environment the NS can set up a pulsar-like 


magnetosphere. Via dipole spin-down, the NS starts power¬ 
ing a highly relativistic, Poynting-flux dominated outflow of 
charged particles (mainly electrons and positrons; see Sec¬ 
tion 4.2.1) or ‘pulsar wind’ at the expense of rotational en¬ 
ergy. This occurs at a time t = fpui.in and marks the beginning 
of Phase II. 

The pulsar wind inflates a PWN behind the less rapidly ex¬ 
panding ejecta, a plasma of electrons, positrons and photons 
(see Section 4.3.1 for a detailed discussion). As this PWN is 
highly overpressured with respect to the confining ejecta en¬ 
velope, it drives a strong hydrodynamical shock into the fluid, 
which heats up the material upstream of the shock and moves 
radially outward at relativistic speeds, thereby sweeping up all 
the material behind the shock front into a thin shell. During 
this phase the system is composed of a NS (henceforth “pul¬ 
sar” in Phase II and III) surrounded by an essentially baryon- 
free PWN and a layer of confining ejecta material. The prop¬ 
agating shock front separates the ejecta material into an in¬ 
ner shocked part and an outer unshocked part (cf. Figure 1 
and 2). While the shock front is moving outward across the 
ejecta, the unshocked matter layer still emits thermal radia¬ 
tion with increasing luminosity as the optical depth decreases. 
Initially, the expansion of the PWN nebula is highly rela¬ 
tivistic and decelerates to non-relativistic speeds only when 
the shock front encounters high-density material in the outer 
ejecta layers. The total crossing time for the shock front is 
typically Afshock = fshock ,out ^pul,in ^ ^pul ,in? where fshock ,OUt 
denotes the time when the shock reaches the outer surface. At 
this break-out time, a short burst-type non-thermal EM signal 
could be emitted that encodes the signature of particle accel¬ 
eration at the shock front. 

Phase III starts at f = fshock,out- At this time, the entire ejecta 
material has been swept up into a thin shell of thickness Agj 
(which we assume to be constant during the following evo¬ 
lution) that moves outward with speed Uej (cf. Figure 2). In 
general, this speed is higher than the expansion speed of the 
baryon-loaded wind in Phase I (Uej.in), as during shock prop¬ 
agation kinetic energy is deposited into the shocked ejecta. 
Rotational energy is extracted from the pulsar via dipole spin- 
down and it is reprocessed in the PWN via various radiative 
processes in analogy to pair plasmas in compact sources, such 
as active galactic nuclei (see Section 4.3.1 for a detailed dis¬ 
cussion). Radiation escaping from the PWN ionizes the ejecta 
material, which thermalizes the radiation due to the optical 
depth still being very high. Only at much later times the ejecta 
layer eventually becomes transparent to radiation from the 
nebula, which gives rise to a transition from predominantly 
thermal to non-thermal emission spectra. We note that for 
reasons discussed in Section 5.6, the total luminosity of the 
system shows the characteristic cx behavior for dipole 
spin-down at late times t ^ tsd, where t^d is the spin-down 
timescale. However, when restricted to individual frequency 
bands, the late time behavior of the luminosity can signifi¬ 
cantly differ from a oc power law. 

As the NS is most likely not indefinitely stable against grav¬ 
itational collapse, it might collapse at any time during the evo¬ 
lution outlined above (see Section 4.4). If the NS is supramas- 
sive, the collapse is expected to occur within timescales of 
the order of ~ t^d, for the spin-down timescale represents the 
time needed to remove a significant fraction of the rotational 
energy from the NS and thus of its rotational support against 
collapse. For typical parameters, the collapse occurs in Phase 
III. However, if the NS is hypermassive at birth and does not 
migrate to a supramassive configuration thereafter, it is ex- 
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Figure 2. Left\ Phase II of the evolution. The high pressure of the pulsar wind nebula (‘nebula’) powered by the spinning-down NS drives a strong shock through 
the ejecta with speed thereby compressing, heating and accelerating the material. Right: Phase III of the evolution. The ejecta material of mass Mej has 
been entirely swept up by the shock into a thin layer of thickness Aej moving outward at speed while the spinning-down NS keeps injecting energy into the 
optically thick pulsar wind nebula. 


peeled to collapse already in Phase I (i.e., on the timescale for 
removal of differential rotation). Moreover, even a supramas- 
sive NS can collapse on timescales shorter than ^ fsd. if the 
collapse is induced by a (magneto-)hydrodynamic instability. 

Assuming that a SGRB and its afterglows are produced by 
a BNS system merging into a long-lived NS, the prompt burst 
can be associated either with the merger itself or with the de¬ 
layed collapse of the remnant NS, as in the recently proposed 
‘time-reversal’ scenario (Ciolfi & Siegel 2015b,a; see Rez- 
zolla & Kumar 2015 for an alternative proposal). Both sce¬ 
narios can be accommodated by the framework of our model. 
In the former case, the lightcurves and spectra we predict after 
the time of merger represent a model for the observed after¬ 
glows. In the latter case, only the emission that follows the 
collapse should be directly compared with the observed after¬ 
glows. Predictions for both scenarios are discussed in detail 
in the companion paper (Paper II). 


3. EVOLUTION EQUATIONS 

According to the three main evolution phases discussed in 
Section 2, our model consists of three sets of coupled ordinary 
differential equations (ODEs), which we think capture the 
main physical ingredients to describe the dynamical evolution 
and EM emission of the post-merger system. These ODEs 
assume spherical symmetry and they are formulated in terms 
of main evolution quantities, such as the extent of certain 
dynamical structures (e.g., the bayon-loaded wind, the PWN 
and the ejecta shell) and their associated energy budgets. In 
summary, the evolution equations we consider are as follows 
(see Table 1 and Eigure 2 for definitions of variables): 

Phase I: baryon-loaded wind 


dJ^ej 

dt 

dSth 




LEuit) + 


di?th,NS 

dt 
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( 2 ) 


Phase II: pulsar wind shock 
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Phase III: pulsar wind nehula 
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(13) 
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Equations (1) and (2) are supplemented with a model to 
describe the baryonic wind emitted from the NS during this 


dt 


-^rad(t) 
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Table 1 

Definition of lab-frame quantities. 


Quantity 

Description 

t 

time 

Rej 

outer radius of the ejected matter 

Rn 

radius of the PWN (Phase II, III) 

Rsh 

radial location of the pulsar wind shock front (Phase II) 

i)„(r, t) 

velocity prohle of the baryon-loaded wind (Phase I, II; 
cf. Section 4.1.1, Equation (30)) 

'^ej,in 

initial expansion speed of the baryonic ejecta material 
(cf. Section 4.1.1, Equation ( 1 9)) 


velocity of the pulsar wind shock front (Phase II; cf. Equa¬ 
tion (58)) 

"llej 

velocity of the ejecta shell in Phase III (cf. Equation (95)) 

Q-ej 

acceleration of the ejecta shell in Phase III (cf. Equation (95)) 

Ash 

radial thickness of the shocked ejecta shell in Phase II 
(cf. Section 4.2.3) 

Aej 

radial thickness of the ejecta shell in Phase III (cf. Sec¬ 
tion 4.3) 

/ej 

function to smoothly ‘switch on/off’ terms as the ejecta ma¬ 
terial becomes optically thin (cf. Equation (102)) 

Rth 

internal energy stored in the ejecta material, available to be 
emitted as thermal radiation 

^th,sh 

internal energy stored in the shocked ejecta matter (Phase II) 

^th,ush 

internal energy stored in the unshocked ejecta matter (Phase 

II) 

internal energy of the PWN 

-S'nth 

Eb 

magnetic energy of the PWN 

dL’th,Ns/dt 

internal energy injected from the NS into the baiyon-loaded 
wind per unit time (Phase I; cf. Equation (33)) 

dRsh/dt 

energy deposited in the shocked ejecta material by shock 
heating per unit time (cf. Equation (63)) 


total internal energy of unshocked ejecta in the volume swept 
up by the shock front per unit time (cf. Equation (62)) 

dL’pwN/dt 

total energy emitted by the PWN per unit time (Phase II, III; 
cf. Equations (57) and (92)) 

i-EM 

EM energy deposited in the baryon-loaded wind per unit time 
(Phase I; cf. Equation (32)) 

-^rad 

luminosity of thermal radiation from the outer surface of the 
ejected material (cf. Equations (34), (74), and (100)) 

-^rad,in 

luminosity of thermal radiation from the ejected material ra¬ 
diated toward the interior (Phase II,III; cf. Equations (73) and 
(101)) 

rad,pul 

luminosity of thermal radiation from the NS surface (Phase 
II,III; cf. Equation (51)) 

^sd 

spin-down luminosity of the pulsar (Phase II, III; cf. Equa¬ 
tion (48)) 


fraction of the total pulsar wind power injected as magnetic 
energy per unit time into the PWN (cf. Sections 4.2.2 and 
4.3.1) 

'Tts 

efficiency of converting pulsar wind power into random ki¬ 
netic energy of accelerated particles in the PWN (cf. Equa¬ 
tions (9), (76) and Sections 4.2.2, 4.3.1) 


phase (see Section 4.1.1). Furthermore, Equations (11)-(15) 
are supplemented with Equations (78) and (79), which model 
the radiative processes inside the PWN and which need to be 
solved at every time step of the evolution equations to deter¬ 
mine the source term di?pwN/df and the emission spectrum 
of the nebula. 

The above evolution equations are formulated in the rest 
frame of the merger remnant (henceforth the “lab frame”; 
see also Appendix B). The various terms appearing in Equa¬ 
tions (1)-(15) are motivated and discussed in detail in Sec¬ 
tion 4. 

4. INGREDIENTS 

4.1. Phase I: baryon-loaded wind 

Phase I starts with a differentially rotating NS that generates 
a baryon-loaded wind, either due to a very strong magnetiza¬ 
tion of the stellar interior (Siegel et al. 2014) and/or as the re¬ 
sult of neutrino emission from the stellar interior (e.g., Dessart 
et al. 2009). This wind can be treated as roughly spherically 


symmetric at distances r > Rmin = 30 km (cf., e.g., Siegel 
et al. 2014; Siegel & Ciolfi 2015c), which we define as our 
inner spatial boundary of the evolution model. 

Such a magnetically and/or neutrino-induced wind is likely 
to dominate baryon pollution around the newly-formed NS on 
the timescales of interest. Dynamical ejecta originating from 
the tidal tails during the merger process will be ejected mostly 
into the equatorial plane and move away from the merger site 
with high (mildly relativistic) velocities (e.g., Davies et al. 
1994; Rosswog et al. 2013; Oechslin et al. 2007; Hotokezaka 
et al. 2013; Bauswein et al. 2013). They are thought to 
undergo r-process nucleosynthesis and to possibly power a 
macronova (aka a kilonova; e.g., Li & Paczyhski 1998; Met¬ 
zger et al. 2010; Barnes & Kasen 2013; Piran et al. 2013; 
Tanaka & Hotokezaka 2013). The more isotropic dynami¬ 
cal ejection originating from shock heating at the contact in¬ 
terface during collision of the NSs and aided by radial oscil¬ 
lations of the double-core structure of the newly-formed NS 
immediately after merger (Hotokezaka et al. 2013; Bauswein 
et al. 2013; Kastaun & Galeazzi 2015) only lasts for a few mil¬ 
liseconds and is thus unlikely to dominate baryon pollution on 
the much longer timescales relevant here. Such a component 
can, however, be accounted for in our wind model by, e.g., 
tuning the mass ejection rate M (see below). 


4.1.1. Wind model 

In order to avoid performing expensive hydrodynamical 
simulations, we describe the generation and evolution of 
the baryon-loaded wind by a simple one-dimensional model. 
In particular, we assume that there are no hydrodynamical 
shocks being formed in the wind. The assumption of spherical 
symmetry is motivated by recent three-dimensional magne¬ 
tohydrodynamic simulations in general relativity, which have 
shown that for the most realistic magnetic field configurations, 
such a wind will be highly isotropic (Siegel et al. 2014; Siegel 
& Ciolfi 2015c). Our simple approach is furthermore justified 
by the fact that the exact details of modeling this wind only 
have a minor influence on the final predictions of our model, 
such as the lightcurves, as the pulsar nebula shock eventually 
sweeps up all the material into a thin shell in Phase II. 

It is conceivable to assume that a fluid element of the wind 
ejected from the NS at time t moves outward with a constant 
velocity v{r,t) = v{t), once it has climbed up the gravita¬ 
tional potential of the NS after a few tens of kilometers. Mass 
conservation requires v{r,t) = M{t)/p{r,t), where p 
denotes the rest-mass density and M{t) the mass injection 
rate into the wind. Hence, in order to satisfy v{r, t) = v{t), 
p{r, t) = A{t) jr^ for matter ejected at time t. This assump¬ 
tion is also supported by recent general-relativistic magneto¬ 
hydrodynamic simulations of magnetically driven winds that 
yielded a p oc density profile and constant ejection veloc¬ 
ities for a constant M over the timescale of tens of millisec¬ 
onds (Siegel et al. 2014; Siegel & Ciolfi 2015c). 

While differential rotation is being removed on a timescale 
/dr, i-e., the gradient of the angular frequency decays as 

|Vfl| oc exp(-//fdr), (16) 

we also expect that the mass ejection rate, therefore the den¬ 
sity profile, and the ejection speed decay on roughly similar 
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timescales; 




time t. 

M{f)=Mm exp (-(TMf/f dr), 

(17) 

m^{rA)= / M 

A{f) = Ain exp(-(Tpf/fdr), 

(18) 

Ji(r) 

v{t)=V^ym exp(-cr„f/fdi), 

(19) 

II 

1 

tJn are of the order of one. This set of assump- 

O'M 


where ctm , cFp 

tions satisfies the ‘no-shocks-requirement’ and immediately 
yields for the outer radius of the wind; 


-^ej(f) — f^min 'Uej,inf- 


( 20 ) 


M{t') dt' 


(27) 


exp -o-M- 


fdr 


-expi -<7m — 
tdr 


where t{r, t) denotes the time a fluid element at r and t was 
sent out from the inner boundary at r = i-S-, t is deter¬ 
mined hy r = — F) + Rmin, or, equivalently, 


Equations (17)-(19) only apply for a magnetically driven 
wind (Siegel et al. 2014). However, typical values for fjr ^ 
1 s agree well with typical neutrino cooling timescales, such 
that a neutrino-driven wind can be additionally incorporated 
in the model by tuning M^- We note that even for mod¬ 
erate magnetic field strengths and realistic rotation periods 
of the NS the magnetically driven wind can easily exceed 
mass loss rates of a few in 10“^ Mq s“^ and thus dominate 
over the other ejection mechanisms. If the neutrino cooling 
timescale is very different from the situation could be dif¬ 
ferent. However, as the pulsar wind shock eventually sweeps 
up all the ejecta material into a thin shell, the following evo¬ 
lution will mostly depend on the total amount of emitted ma¬ 
terial Mq (X Minfdr- Hence, by varying fdr and Min the above 
model for mass ejection can effectively also accommodate a 
dominant neutrino-induced wind. 

Simple scaling arguments can fix two of the free parame¬ 
ters aM,o'p,(7v that control the decay timescales relative to 
fdr- Magnetic winding in the stellar interior converts rota¬ 
tional energy into magnetic energy (mostly into toroidal field 
strength) at a rate cx | VHp (e.g., Duez et al. 2006; Siegel et al. 
2013). This magnetic energy is then available to be dissipated 
into kinetic and EM energy of a wind driven by the built-up 
of magnetic pressure in the stellar interior (Siegel et al. 2014). 
Consequently, 


{t — t) exp 



r — R 


min 


^ej,in 


= 0 . 


(28) 


It is straightforward to show that Equation (28) has exactly 
one root in [0, t], such that t{r, t) is well defined. At any time 
t the density profile of the wind is then given by 


1 dm„{r,t) 


The corresponding velocity profile is 






exp —a. 


’ fdr 


(29) 


(30) 


and the total amount of ejected mass at time t is given by 
Mej(f)=mw(i?ej(f),f) 


= —M„ 

O'M 


1 - exp( -(Tm^ 
Cdr 


(31) 


4.1.2. Properties of the ejected matter 

Injection of electromagnetic energy — The wind will be en¬ 
dowed with a Poynting flux of luminosity (see Siegel et al. 
2014 and Section 4.1.1) 


TEM(f) — 10 


,48 


B 


VlOiSQ 


R. y 

10® cm ) 


Pc 

10-4 S 


-1 


-Mv"^ + Lem oc |VH|" 


( 21 ) 


X exp 



ergs , 


(32) 


where Lem denotes the EM luminosity corresponding to the 
Poynting flux carried by the wind (cf. also Section 4.1.2). 
Prom Equations (16) and (21) we conclude that Lem fx 
exp(—tTif/idr), with ctl = 2, and 


(Jm T 2(j^ — 2. 


( 22 ) 


At r = Rmin, we also need to satisfy v{t) = 
M{t)/Altp[Rmin,t)Rmin- U^ing Equations (17)-(19), this 
yields 


<J M — CTp -f 


(23) 


where cti, = 2. While fdr is a free parameter of our model, the 
total magnetic field strength in the outer layers of the NS, B, 
the equatorial radius i?e, as well as the central spin period 
of the differentially rotating NS can be extracted from numer¬ 
ical relativity simulations. The wind itself is highly turbulent 
in the vicinity of the NS (Siegel et al. 2014) and magnetic dis¬ 
sipation will therefore be very effective. We assume that this 
Poynting flux is dissipated and thermalized in the ejecta mat¬ 
ter (due to the very high optical depth) on the timescales of 
interest, such that this energy is trapped and thus appears as a 
source term in Equation (2). 


and Ain = Mm/A-KVeyin- According to Equations (22) and 
(23), only one a-parameter can be chosen independently. The 
additional requirement of ctm, (j„ all being non-negative 
limits these parameters to the following ranges; 

ctmG[1,2] (24) 

tTpG [0.5,2] (25) 

(t^G [0,0.5]. (26) 

In order to compute the density profile p(r, f) of the wind 
resulting from the mass ejection model (17)-(19), we first 
compute the total mass contained in a volume of radius r at 


Injection of thermal energy — The material ejected from the 
NS surface as probed by numerical relativity simulations typ¬ 
ically has a very high specific internal energy Cej.NS.in that cor¬ 
responds to a temperature Tq^NS.in of the order of tens of MeV 
(e.g., Siegel et al. 2014; Kastaun & Galeazzi 2015). We there¬ 
fore source the internal energy of the ejected material in Equa¬ 
tion (2) by the corresponding injection rate 


dLth.NS 

dt 


^ej,NS,in.Min CXp 




(33) 


where we have assumed that the NS matter cools down due to 
neutrino cooling on a timescale 
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Thermal radiation from the wind — Due to the very high optical 
depth of the wind material, the radiation originating from the 
surface of the expanding wind will be predominantly thermal 
in Phase I. This leads to radiative losses in Equation (2) with 
a luminosity 

L,,d(t)=47ri?2(f)aTjr(f), (34) 


where a denotes the Stefan-Boltzmann constant. The effec¬ 
tive temperature of the wind, Teff, is given by 




16 T^it) 

3 AT(f) + 1 


(35) 


Here, 


pRa] {i) 

Ar(f) = K / /9w,t(r)dr 

jRUt) 


(36) 


is the optical depth of the wind with associated diffusion 
timescale 

tMt) = + 1], (37) 

c 

Adding unity here and in the definition of Teff takes a possible 
transition to the optically thin regime into account, k is the 
opacity of the ejecta material and Rin{t) denotes the effective 
inner radius of the ejected material, defined by the condition 


K 





= 1 . 


(38) 


We describe the wind material as a mixture of ideal gas with 
adiabatic index Tej = 4/3 and radiation. The temperature 
T(t) of the wind in Equation (35) at time t is thus obtained as 
the root of the following equation: 


T{tf + 


1 fcfi 

Tej - 1 TOpO 




Eth{t) 


(39) 


where is the Boltzmann constant, rrip denotes the proton 
mass, a is the radiation constant, T4j(f) = (4/3)7r[i?4(f) — 
the volume of the ejecta material, and p„{t) = 
— mw(i?in(t), f)]/14j(t) its mean density. 


4.2. Phase II: pulsar wind shock 
4.2.1. Pulsar properties 

‘Ignition ’— As differential rotation is being removed the 
NS settles down to uniform rotation and mass ejection is 
suppressed according to the wind model discussed in Sec¬ 
tion 4.1.1. This eventually creates the conditions to build up a 
pulsar magnetosphere. As a criterion to ‘switch on’ the pulsar, 
we employ p(i?niin,i) = (Ain/ii^jJ exp(-crpf/tdr) < Pent, 
which translates into a time 

= (40) 

p \ Afn J 

Eor the critical density we typically choose the rest-mass 
density for Iron ions corresponding to the Goldreich-Julian 
charge density, peUt Petpe] ~ 6 x 10“®gcm“^, assuming 
typical magnetic fields of 10^® G and rotational periods of the 
pulsar of ~ 1ms (Goldreich & Julian 1969). We note, how¬ 
ever, that as the density in the vicinity of the NS decreases 
exponentially, varying the choice of pciit even by orders of 
magnitude does not influence the model evolution noticeably, 
as long as the value is sufficiently small. 


Initial spin and magnetic field — Once the NS has settled down 
to uniform rotation, it is expected to maintain roughly the 
same level of magnetization on the timescales of interest. In 
Phase I, strong toroidal fields have been built up through mag¬ 
netic winding and possibly instabilities like the magnetorota- 
tional instability (Duez et al. 2006; Siegel et al. 2013; Siegel 
& Ciolfl 2015b; Kiuchi et al. 2014), such that the strength of 
the dipolar component of the poloidal held at the pole Bp is 
only a fraction rjB^ of the total magnetic held strength in the 
outer layers of the NS, 

Bp = Vb^B. (41) 

This dipole component is the relevant quantity for dipole spin- 
down emission (see next paragraph). Knowing the initial ro¬ 
tational energy Trot.NS.in of the differentially rotating NS from 
numerical relativity simulations, we can infer the initial spin 
period of the pulsar Ppui.in at fpuUn from its initial rotational 
energy Trot.pui.in: 

Prot.pul.in — Prot,NS,in BeM Pkin.w Prot.ej • (42) 


Here, 


Eem — 


LEuit) dt 


(43) 


is the rotational energy dissipated into Poynting flux during 
Phase I, 


^ (Ipul.in) 

Pkin.w = ;,47r / Pw.v.nW^'w(AVun)»'^dr (44) 

^ I Rmia 


is the rotational energy dissipated into kinetic energy of the 
radially expanding wind in Phase I, and 


E, 


'rot,ej 


-^ej (^pul,in) 

AfNS.in 


E, 


rot,NS,in 


(45) 


approximates the (subdominant) amount of initial rotational 
energy carried away by the ejected matter. The initial mass 
Afxs.in of the NS at birth is known from numerical relativity 
simulations. The initial spin period is then given by 

_ 1 

Ppul,rn = 27r , (46) 

where /pui is the moment of inertia of the pulsar. 

Spin-down luminosity — Due to unipolar induction charged 
particles are continuously extracted from the pulsar surface 
(Goldreich & Julian 1969). This primary particle current ini¬ 
tiates copious pair production, which populates the magne¬ 
tosphere with a nearly force-free plasma and drives a highly 
relativistic, magnetized (Poynting-flux dominated) outflow of 
particles (mostly electrons and positrons), referred to as a pul¬ 
sar wind. The associated Poynting-flux luminosity is given by 
(Spitkovsky 2006; Philippov et al. 2015) 

/ry \4 D2 d6 

Lsd,m = p 4 {ko + h sin^ x), (47) 

where kg — 1.0 ± 0.1, ki = 1.1 ± 0.1, x the inclination 
angle of the dipole component of the magnetic held with re¬ 
spect to the rotation axis of the pulsar, and i?pui — Re denotes 
the radius of the pulsar. This Poynting flux is extracted from 
the pulsar via the magnetic held at the expense of rotational 
energy Trot.pui- Observations of well studied objects like the 
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Crab pulsar suggest that only ~ 1% of the observed spin- 
down energy is directly radiated away by the pulsar (Buhler 
& Blandford 2014), such that essentially the entire spin-down 
energy is carried by the pulsar wind leaving the light cylinder. 
Hence, we can set i^rot.pui = which results in slowing 
down the pulsar spin Ppui(f) oc -\- t/tsd and yields a spin- 
down luminosity of 


where 


-^sd(^) — -^sd,in | 

" t-tp,Un\ " 

, tsd ) ’ 

(48) 

^sd ~ 

-£'rot,pul,m 

(49) 

-^sd,m 


is the spin-down timescale. In the following, we assume an 
aligned rotator, i.e., x = 0, as a non-zero inclination angle 
would result in a change of Lsd of order unity, which can be 
absorbed, e.g., in the parameter that sets the initial field 
strength Bp (see Equation (41)). 


Metzger et al. 2014; Metzger & Piro 2014): we formulate the 
energetics of the PWN in terms of a balance equation for the 
total internal energy E^th of the nebula (photons and random 
kinetic energy of the particles; cf. Equation (9)). The internal 
energy of the PWN is sourced by the spin-down luminosity 
Lsd (Equation (48)), thermal radiation from the pulsar surface 
Trad,pul (Equation (51)) as discussed above (which can rep¬ 
resent a significant contribution at early times), and thermal 
radiation emitted by the hot ejecta matter toward the interior 
with luminosity Lrad.in (Equation (73)). We assume that only 
a fraction pxs of the spin-down luminosity and the thermal 
radiation from the pulsar surface are dissipated into random 
kinetic energy of the particles in the nebula. This parameter 
reflects our incomplete knowledge about the efficiency with 
which the bulk energy of the pulsar wind is dissipated in the 
nebula (see Section 4.3.1 for more details). 

Rapid expansion of the PWN occurs at the expense of p AV 
work, 

TE/pdV AVji E'nth rli^n 


Emission of thermal energy — Immediately after the BNS 
merger the NS matter is very hot with temperatures of a 
few tens of MeV. It cools down via neutrino emission on a 
timescale < Is. This thermal energy reservoir can be 
tapped via ejection of material in Phase I (cf. Equation (33)) 
and by thermal radiation from the pulsar surface in Phase 
II and III. Thermal radiation from the pulsar can help initi¬ 
ate secondary-particle cascades based on the primary charged 
particles extracted from the pulsar surface. The associated en¬ 
ergy is thus reprocessed in the magnetosphere and will be car¬ 
ried away by the pulsar wind. However, given typical neutrino 
cooling timescales, this thermal energy contribution will be¬ 
come energetically negligible very soon with respect to, e.g., 
the spin-down luminosity. Therefore, we adopt a simple cool¬ 
ing model for the pulsar surface. 


where 


_ Enth , 

“ SK Stt 


(53) 


is the pressure inside the nebula with volume 14 = (4/3)7ri?j^ 
and radius R„. Here we have assumed that the particles inside 
the nebula are relativistic and collisionless (cf. Section 4.3.1), 
such that they form an ideal gas with adiabatic index 4/3. 
The second pressure contribution on the right hand side of 
Equation (53) is caused by the uniform and isotropic magnetic 
field that we assign to the nebula. Its field strength at time 
t is inferred from the total magnetic energy 


Stt 


(54) 


Tm{t) = TNs.inexp^—— ^ , (50) 

with associated thermal radiation of luminosity 

Lmd,pui{t) = 47riip„icrT4(f). (51) 

Here, TNs.in is the initial typical surface temperature of the NS 
corresponding to the specific internal energy Cej.NS.in (cf. Sec¬ 
tion 4.1.2). As the thermal energy is reprocessed in the mag¬ 
netosphere, we add Lrad.pui as a source term to the evolu¬ 
tion equation (9) and as an input for the nebula in Phase III 
(cf. Section 4.3.1). 


of the nebula, which we evolve according to Equation (10). 
The parameter riB„ in Equation (10) controls the level of mag¬ 
netization of the nebula (see also Section 4.3.1). We note that 
the efficiency parameters have to satisfy t/ts + < 1- 

Eurthermore, the nebula loses energy through irradiation of 
the surrounding ejecta matter. The internal energy of the neb¬ 
ula is radiated away on the diffusion timescale through various 
radiative processes that we do not model in detail in Phase II 
(see, however. Section 4.3.1). We may write the associated 
luminosity as Lnth(f) = £^nth(i)/4iff,n(i), where 

tdm,n(t) = ^ [1 -f ATT(f)] (55) 


4.2.2. Simple model for the expanding pulsar wind nebula 

The newly-formed pulsar wind (Section 4.2.1) leaves the 
magnetosphere with relativistic velocities and inflates a PWN 
behind the less rapidly expanding ejecta material (e.g.. Ken¬ 
nel & Coroniti 1984a; see Section 4.3.1 for more details). The 
PWN is highly overpressured with respect to the surrounding 
ejecta matter and thus drives a strong hydrodynamical shock 
into the material (see Section 4.2.3), which, in turn, leads to a 
rapid expansion of the PWN. The exact physical description 
of such a highly dynamical PWN is complex. However, the 
main energetical features that govern the overall dynamics of 
the system (which is what we are interested in here) can be 
captured using a simple approach inspired by recent dynami¬ 
cal PWN models (e.g., Gelfand et al. 2009; Kotera et al. 2013; 


is the photon diffusion timescale of the nebula due to Thom¬ 
son scattering of photons off thermal electrons and positrons. 
Here, 


ArT(f) 


4yo-TTsd(f) 

TTRn{t)m^C^ 


(56) 


is the optical depth to Thomson scattering, with aj the Thom¬ 
son cross section, rrie the electron mass, and Y the pair yield 
of the nebula (Lightman & Zdziarski 1987). Given the very 
high electron compactness parameter /g ^ 1 during Phase II 
(see Section 4.3.1), the pair yield is given by K ~ 0.1 (Svens- 
son 1987; Lightman & Zdziarski 1987). We have added the 
light crossing time in Equation (55) to account for a possible 
transition to the optically thin regime. 
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Ionization of ejecta matter due to nebula radiation can cause 
a non-zero albedo A, which, in turn, can be frequency de¬ 
pendent, such that the radiative losses of the nebula may be 
written as (cf. also Metzger et al. 2014; Metzger & Piro 2014) 

= /beam(t) J “ A{t, i/)]Lnth(f, v) (57) 

with f L„±(t, ly) diy = Lnth(f). Due to relativistic expansion 
of the PWN with velocity Vn = dRa/dt (cf. Equation (5)), 
relativistic beaming will prevent a sizable fraction 1 — /beam 
of photons to escape from the nebula. Therefore, we have ad¬ 
ditionally included a factor /beam(i) = 1 — in Equa¬ 

tion (57). As the radiative losses given by Equation (57) are 
thermalized in the surrounding ejecta shell, the resulting light 
curve neither depends on the exact spectral shape of nebula ra¬ 
diation Tnth(i) v) nor on the spectral shape of A{v^ t). There¬ 
fore, we typically set A{t, v) = A{t) and we need not assume 
a spectral shape for Lnth(f)- 


the shocked ejecta layer is maintained. Eor a mixture of ideal 
gas with Tej =4/3 and radiation, this results in the following 
constraint: 

i?n=i?sh^. (61) 

r^nth 

This conditions defines dAsh/df and couples Equations (5) 
and (9) to make them implicit evolution equations for and 

-E'nth- 


Shock related energetics — The shocked and unshocked ejecta 
material are assigned internal energies E'th.sh and i^th.ush, re¬ 
spectively, which are evolved separately (cf. Equations (6) 
and (7)). As the shock propagates across the ejecta matter in¬ 
ternal energy is transfeiTed from the unshocked to the shocked 
fluid part, which amounts to 


dE/th.vol 


Ti^th.ush 


(7?sh + di?sh)^ 

R%-Rl 


(62) 


4.2.3. Shock dynamics and ejecta properties 

The high nebula pressure (53) drives a strong hydrodynamic 
shock into the ejecta matter. The shock front (denoted by the 
radial coordinate i?sh) divides the ejecta matter into shocked 
ejecta with internal energy 77*,sh and unshocked ejecta with 
internal energy i^th.ush (cf. Eigure 2). The PWN with radius 
i?n is separated from the shocked ejecta material by a contact 
discontinuity. We denote the radial thickness of the shocked 
ejecta layer by Ash = Rah — Rn- 


during a time dt. Eurthermore, as the fluid is swept up by the 
shock it is heated at a rate 


where 


Ae 


d7/sh 

dt 


Ae 


dmsh 

dt 


1 Pn _ 

Tej + 1 Pr ^ 


2 



4rcj Pn 

(rej-rr)2 prc^ 


2 


(63) 

(64) 


Shock propagation — In order to avoid performing hydrody- 
namical simulations, we describe the shock dynamics by a 
simple model (Equations (4) and (5)). According to the 
special-relativistic Rankine-Hugoniot conditions, the shock 
speed as measured in the lab frame is given by (see Ap¬ 
pendix A) 


t^sh(7) = 


tlw(77sh(f)5 7) -f r^sh,R(^) 

1 + Vw(7?sh(f),f)t^sh,R(f)/c^’ 


where 


(58) 


tIsh.R — 


Pn 


Pr 


t — (t -\ _/a_ 

^ PL + rej-r PLCV + 


Pn 

PRC^ 


(59) 


is the shock velocity in the frame comoving with the un¬ 
shocked fluid. Here, pR(f) = Pw,t(7?sh(f)) is the density of 
the unshocked fluid at the shock front and 


Pl 


Tej + 1 
Dej-l 




1+4 


Pn lej 
Prc 2 (Tej + 1)2 



(60) 


is the density of the shocked fluid at the shock front (see Ap¬ 
pendix A). A special-relativistic framework for shock propa¬ 
gation is required here due to the very low baryon densities in 
the surrounding of the pulsar (cf. Equation (40)) and the very 
high pressure of the pulsar nebula (Equation (53)), which re¬ 
sult in highly relativistic initial shock propagation speeds. The 
shock front decelerates as it propagates outward into higher 
density regions further out from the pulsar and eventually 
slows down to non-relativistic or mildly relativistic velocities 
once it reaches the outer surface of the ejecta material. As the 
shock front is evolved using Equation (4), the thickness Ash of 
the shocked ejecta layer and thus i?n (cf. Equation (5)) is ad¬ 
justed such that pressure equilibrium between the nebula and 


is the jump in specific internal energy across the shock front 
(see Appendix A) and 


dwsh 

dt 


47ri?shPR 



[t^sh t^w(77sh5 *)] 


(65) 


is the amount of material swept up by the shock per unit time. 
The latter expression is obtained by noting that the amount of 
material in a spherical volume changes by 


dm = 



dirr'^pv 



df 


( 66 ) 


if the radius of the volume r is altered by dr during a time 
dt, assuming that the fluid is radially moving outward with 
velocity v. 

Eor further reference (see Section 4.3.2), we also note that 
the shock deposits kinetic energy in the shocked ejecta ma¬ 
terial. The jump Ackin in specific kinetic energy across the 
shock front as measured in the lab frame is given by 



where = Vv,{Rsh{t),t) is the velocity of the wind in 

the unshocked fluid part at the shock front as measured in the 
lab frame. Eurthermore, 


Vhit) + Vsh{t) 

1 + vi^{t)v^h(t)/A 


( 68 ) 


is the velocity of the shocked fluid at the shock front as mea¬ 
sured in the lab frame, with vl being the velocity of the 
shocked fluid at the shock front in the frame comoving with 
the shock front (see Appendix A). 
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Thermal emission — As in Phase I, we assume that the ejecta 
matter consists of a mixture of ideal gas and photons. The 
corresponding fluid temperatures Tgh and Tush are given by 


T 


sh,ush 


-f 


r^j 1 TTIp (X 


_ ry-i '£^th,sh,ush ^ flZCW 

Psh,ush-tsh,ush-777- = U, (09) 


aV^h. 


,ush 


where psh(t) = m^{Rsh{t),t)/Vsh{t) and push(t) = 
[mv,{Rej{t),t) — TOw(Tsh(f),f)]/Kish are the mean densities 
of the shocked and unshocked ejecta, respectively, with asso¬ 
ciated volumes 14h = 47r(i?^|, — i?u)/3 and Ksh = 47r(i?4 — 
The effective temperatures are then defined by 


T 


eff,sh,ush 


it) 


16 ^sh,ush(^) 

3 ATsh,ush(f) + 1’ 


with optical depths 


ATsh(f) = KPsh(i)Ash(t) 


and 


ATush(f) = K 



/Ow,t(T) dr. 


(70) 


(71) 


(72) 


These effective temperatures give rise to thermal radiation 
emitted from the inner surface of the shocked ejecta toward 
the interior with luminosity 


7^rad.in(i) = in Rl{t) (73) 

and from the outer surface of the ejected material with lumi¬ 
nosity 

LUt) = 47rT^j(f)aTu7ush(i)- (74) 

While the radiation associated with Trad leaves the system 
(cf. Equation (7)), Trad,in is reabsorbed by the nebula, its en¬ 
ergy is reprocessed and eventually reemitted into the ejecta 
material (cf Equations (6) and (9)). 

Shock break-out — At the time when the shock front reaches 
the outer ejecta surface a short transient signal can be pro¬ 
duced in addition to the thermal radiation from the ejecta. 
This signal can carry the signature of particle acceleration at 
the shock front and it can have a non-thermal spectrum reach¬ 
ing higher maximum energies than the background thermal 
emission; therefore, it could be observable in the hard X- 
ray and gamma-ray bands. Within the time-reversal scenario 
(Ciolfi & Siegel 2015b), this transient signal might provide a 
convincing explanation for the early precursors observed up 
to ~ 10^ s prior to some SGRBs (Troja et al. 2010). 

4.3. Phase III: pulsar wind nebula 

Phase III starts at fshock.out, when the shock front reaches 
the outer surface of the ejecta material, i?sh(fshock,out) = 
7?ej(Thock,out)- At this time the shock has swept up all the 
ejecta material of mass Mej = Mej(fpui,in) into a thin shell 
of thickness Aq = Ash(fshock,out) (cf. Eigure 2), which we 
assume to be constant during the following evolution. There¬ 
fore, we set di?n/dt = di?ej/df for the evolution of the PWN 
radius R^ in Phase III (cf. Equation (13)). The following two 
subsections detail the physical description of the PWN and 
the ejecta shell in this phase. 


4.3.1. Physics of the pulsar wind nebula 

Basic picture and assumptions — PWNe are traditionally stud¬ 
ied in the context of supernovae-bom pulsars (see Gaensler & 
Slane 2006, Kargaltsev et al. 2015, and Biihler & Blandford 
2014 for reviews). Only very recently, PWNe have also been 
considered in the BNS merger scenario (Bucciantini et al. 
2012; Metzger & Piro 2014). We study such a system here 
motivated by the fact that long-lived NSs turn out to be very 
likely outcomes of BNS mergers (cf. Section 1) and by the 
fact that the baryonic ejecta (cf. Section 4.1) naturally pro¬ 
vide the conditions for such a nebula to form (such as pro¬ 
viding the necessary confining envelope). Although the exact 
physical conditions for a PWN in a BNS merger scenario are 
different from those in a supernova event, it is conceivable 
that qualitatively similar morphologies and dynamics emerge. 
Eor Phase III, we develop a physical description for a PWN 
in the context of BNS mergers inspired by studies of super¬ 
novae and radiative processes in active galactic nuclei. While 
spectra of ordinary PWNe from radio to gamma-ray energies 
can be explained in terms of synchrotron and inverse Comp¬ 
ton emission (for the best studied case, the Crab nebula, see, 
e.g.. Kennel & Coroniti 1984b; Atoyan & Aharonian 1996; 
Volpi et al. 2008; Olmi et al. 2014; Porth et al. 2014), the BNS 
merger scenario requires a more detailed treatment of radia¬ 
tive processes due to the presence of the strong photon field 
from the confining hot ejecta material. In order to achieve 
this and in order to avoid performing expensive MHD simu¬ 
lations, we neglect spatial variations and adopt a ‘one-zone’ 
model for the PWN, in which the nebula is considered to be 
homogeneous and isotropic. The resulting description of a 
PWN we propose for Phase III is much more detailed than in 
Phase II for two reasons. First, Phase III is less dynamical and 
a quasi-stationary description of the PWN can be adopted that 
facilitates a detailed treatment of radiative processes. Second, 
as the ejecta material surrounding the PWN becomes trans¬ 
parent to the nebula radiation at some point, radiation from 
the nebula itself is not reprocessed and thermalized anymore 
and its spectrum becomes directly observable. We therefore 
employ a formalism that predicts the nebula spectrum in a 
self-consistent way. 

The newly formed pulsar loses rotational energy at a rate 
Tsd given by Equation (48) through a highly relativistic mag¬ 
netized particle outflow, referred to as a pulsar wind (cf. Sec¬ 
tion 4.2.1). As in Phase II (cf. Section 4.2.2), we also assume 
that thermal radiation from the pulsar surface with luminos¬ 
ity Tjad.pui (cf- Equation (51)) is deposited in this pulsar wind; 
this can represent a significant energy contribution at early 
times when the NS is still very hot. The pulsar wind inflates 
a PWN behind the less rapidly expanding ejecta material, a 
bubble of radiation and charged particles (mainly electrons 
and positrons), which is separated from the pulsar wind by a 
termination shock at a distance (Gaensler & Slane 2006) 


Rts = 



(75) 


from the pulsar, where the ram pressure of the pulsar wind 
equals the intrinsic pressure pn of the nebula; f denotes the 
fraction of a sphere covered by the wind. 

Magnetosphere models predict that the pulsar wind leaving 
the magnetosphere at the light cylinder i?LC = cPpui/27r is 
highly magnetized, ^ i (e.g., Arons 2012). Here, 
denotes the magnetization parameter, defined as the ratio of 
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Poynting flux to particle energy flux (e.g.. Kennel & Coro- 
niti 1984a). One-dimensional and two-dimensional axisym- 
metric MHD models only reproduce the observed properties 
of the Crab nebula well if cTm ^ 1 just upstream of the ter¬ 
mination shock (e.g., Kennel & Coroniti 1984a,b; Atoyan & 
Aharonian 1996; Volpi et al. 2008; Olmi et al. 2014). Dis¬ 
sipating the Poynting flux into kinetic energy of the parti¬ 
cle flow between the light cylinder and the wind termination 
shock is known as the tr-problem. Recent three-dimensional 
MHD simulations, however, suggest a solution the cr-problem 
(Mizuno et al. 2011; Porth et al. 2014). The simulations by 
Porth et al. (2014) are able to reproduce the morphology of 
the Crab nebula even for stronger magnetizations of up to 
(Jn, ~ few, thanks to magnetic dissipation inside the PWN. 
Enhanced turbulence downstream of the termination shock 
in the 3D simulations efficiently dissipates magnetic energy, 
such that the magnetic energy of the nebula is only a very 
small fraction of the internal energy despite a high magneti¬ 
zation of the pulsar wind at the termination shock. This is in 
agreement with previous observations and spectral modeling 
that infer the magnetic energy of the Crab nebula to be less 
than « 3% of the total energy budget (Biihler & Blandford 
2014, Porth et al. 2014, and references therein). In our model 
of the PWN, we control the level of magnetization of the neb¬ 
ula by introducing a factor r]B„, which specifies the amount of 
inflowing magnetic energy per time in terms of the total power 
of the pulsar wind: Eb = + irad,pui(f)] (cf. Equa¬ 

tions (10) and (15)). The magnetic energy Eb then defines 
an average magnetic field strength of the nebula through 
Equation (54). 

At the wind termination shock the wind plasma is deceler¬ 
ated and heated, and efficient conversion of flow energy into 
particle acceleration is thought to take place. The acceler¬ 
ated particles, which are characterized by non-thermal energy 
spectra, are then advected with the flow while cooling down 
due to synchrotron emission. Observations of the Crab nebula 
synchrotron emission indicate that the conversion of pulsar 
wind energy into accelerated particles is > 10% (Kennel & 
Coroniti 1984a; Biihler & Blandford 2014; Olmi et al. 2015; 
we henceforth denote this efficiency by t/ts)- While diffusive 
shock acceleration at the termination shock is now thought 
to be rather unlikely to occur in the Crab nebula and other 
PWNe (Arons 2012; Buhler & Blandford 2014), enhanced 
turbulence and magnetic dissipation downstream of the ter¬ 
mination shock (as found in recent 3D MHD simulations) 
might provide the required non-thermal particle acceleration 
site. While the processes of particle acceleration still remain 
unclear, assuming that the pulsar wind consists of electrons 
and positrons, and that these particles are reaccelerated into 
a power-law spectrum at the termination shock, ID, 2D, and 
3D MHD models of the Crab nebula have been able to repro¬ 
duce the observed photon spectra very well (Kennel & Coro¬ 
niti 1984b; Atoyan & Aharonian 1996; Volpi et al. 2008; Ca¬ 
mus et al. 2009; Olmi et al. 2014; Porth et al. 2014; Olmi et al. 
2015). 

Eollowing this approach we consider a pulsar wind consist¬ 
ing of electrons and positrons and assume that these particles 
are injected into the nebula with a power-law spectrum, al¬ 
though more complex injection spectra can easily be accom¬ 
modated by our model. In view of our ‘one-zone’ descrip¬ 
tion of the nebula, we assume that electrons and positrons 
are continuously injected uniformly throughout the nebula, 
instead of being injected at the termination shock and then 


being advected with the flow. We specify the injection of par¬ 
ticles by the following dimensionless compactness parameter 
(in analogy to Guilbert et al. 1983; Svensson 1987; Lightman 
& Zdziarski 1987): 


h{t) — ^ ^Ts[i'sd(i) + i'rad,piil(f)] 


iira-YRl (^) 
3c 



l)d7. 


(76) 


where = Qo{t)j is the number of particles in¬ 

jected per unit time per unit normalized energy 7 = e/rrieC^, 
where e denotes the particle energy, and per unit volume of 
the nebula. The parameter ? 7 ts defines the aforementioned ef¬ 
ficiency of converting pulsar wind power into random kinetic 
energy of accelerated particles. The power-law injection pa¬ 
rameters r]js, 7 max, and Te are model input parameters, which, 
in the case of observed PWNe such as the Crab, are usually 
determined by comparing simulated emission with observa¬ 
tional data. 

In contrast to ordinary PWNe, where intrinsic photon 
sources inherent to the pulsar-nebula system are absent and 
only background photons from the cosmic microwave back¬ 
ground (CMB) and potentially from local dust and starlight 
are ‘injected’ into the nebula, we need to take additional pho¬ 
ton sources into account in our description of the PWN that 
generate a strong photon field. The hot shock-heated ejecta 
matter confining the PWN injects thermal photons with lumi¬ 
nosity Lrad.in (cf. Section 4.3.2, Equation (101)). For com¬ 
pleteness, we also include a thermal input spectrum Lcmb 
from the CMB (see below), once the ejecta shell has become 
optically thin. We specify the resulting photon injection into 
the nebula in terms of the following dimensionless compact¬ 
ness parameter (in analogy to Lightman & Zdziarski 1987): 


lph{t) = 


(7j 


meC^Rn{t) 

ATTajRlit) 


3c 


[7vrad,in(i) + /ej (f )7vCMB (f)] 

/ ho{x,t)xdx, 


(77) 


where no denotes the combined number of photons injected 
per unit time per unit dimensionless energy x = hn/m^c^ 
per unit volume of the PWN by the aforementioned sources 
(cf. Equation (80)). Here, h denotes the Planck constant and 
1 / is the photon frequency. Photon energies are assumed to 
range between x^m and x^ax, which we define in Section 5.2. 
The function /ej is designed to ‘switch on’ the (subdominant) 
CMB contribution once the ejecta material becomes optically 
thin (cf. Equation (102)). 

Radiative processes in the PWN reprocess the injected par¬ 
ticles and photons and determine the radiative losses of the 
nebula that are transferred to the surrounding ejecta (denoted 
by di?pwN/df in Equation (14)) and the associated emergent 
spectrum. In particular, they determine how spin-down en¬ 
ergy of the pulsar is transmitted to the ejecta material to be 
radiated away from the system eventually. In contrast to or¬ 
dinary PWNe, the photon field in our case is typically very 
strong, such that, expressed in terms of the compactness pa¬ 
rameters defined above, Zph can become comparable to or 
even larger. The photon and particle spectra inside the neb¬ 
ula thus become highly coupled and the computations become 
intrinsically non-linear. In order to determine those spectra 
in a self-consistent way under the combined effects of syn¬ 
chrotron losses, (inverse) Compton scattering off thermal and 
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non-thermal particles, pair production and annihilation, and 
photon escape, we employ, extend and modify a formalism 
developed for pair plasmas in compact sources such as active 
galactic nuclei (Lightman & Zdziarski 1987). 

Balance equations — Given the setup described above, the 
physical processes determining the PWN properties and emer¬ 
gent spectra are remarkably similar to those in pair plasmas of 
active galactic nuclei. Theoretical models to compute detailed 
emergent spectra for active galactic nuclei have been devel¬ 
oped by many authors in the past (e.g., Lightman & Zdziarski 
1987; Coppi 1992; Belmont et al. 2008). Consistent with our 
assumptions, we adopt the approach of Lightman & Zdziarski 
(1987) (henceforth LZ87), but extend and modify it to, e.g., 
include synchrotron losses. We shall briefly outline the for¬ 
mulation of the equations here with emphasis on the modifi¬ 
cations and refer the reader to LZ87 for more details on all 
other aspects. 

According to our ‘one-zone’ model we assume that parti¬ 
cles and photons are uniformly injected throughout the spher¬ 
ical volume of the PWN with radius R„ and associated com¬ 
pactness parameters 4 and Zph (Equations (76) and (77)). At 
a time t and assuming quasi-stationarity, the number densi¬ 
ties per unit energy of photons, n{x), and non-thermal par¬ 
ticles, N{j), are determined by the following set of highly 
non-linear, coupled integro-differential equations; 

0 = flQ -f TIq flsyn 

-^n(ATc'^ -f At^j) - fiesc, (78) 

Jtn 

0 = g(7)+-P(7) + ^c,syn(7)- (79) 

Photons at energy x are produced at rates ho, hp,, 
nj, and fisyn via injection (cf. Equation (80)), e^-pair anni¬ 
hilation, Compton scattering off non-thermal electrons and 
positrons (cf. Equation (9) in LZ87), Compton scattering 
off thermal electrons (cf. Equations (22) and (23) in LZ87; 
see below), and synchrotron cooling of particles (cf. Equa¬ 
tion (87)), respectively. Eor the pair annihilation term hp, 
we use Equations (7), (8), (11), (15)-(18) of Svensson (1983). 
Photons at energy x are lost due to Compton scattering off 
non-thermal particles with optical depth (cf. Equa¬ 

tion (10) in LZ87), due to e^-pair creation with optical depth 
Ar^y^ (cf Equation (11) in LZ87), and by escaping from the 
nebula at a rate nesc (see Equation (21) of LZ87). Addition¬ 
ally, photons at energy x are lost via Compton scattering off 
thermal particles (see below), which is accounted for in nj. 
Photons not ‘absorbed’ by one of these processes are still 
impeded from escaping the nebula by Thomson scattering 
off thermal particles (see below) with scattering depth Atj 
( cf. Equation (20) in LZ87). 

The rate of change of particles at Lorentz factor 7 is given 
by injection into the nebula Q (cf. Equation (76)), e^-pair cre¬ 
ation P (cf. Equation (13) in LZ87), and non-thermal Comp¬ 
ton scattering and synchrotron losses, which are denoted by 
Ncpyn (see below). Particles do not escape from the neb¬ 
ula, i.e., they are assumed to be trapped, e.g., by the mag¬ 
netic held Bn of the nebula. Once cooled down to 7 ^ 1 
and before annihilating, particles are assumed to thermal- 
ize and to form a distinct thermal population described by a 
Maxwell-Boltzmann distribution with dimensionless temper¬ 
ature 0e = koTn/rrinC^ <C 1. The number density of this 
thermal population is determined in terms of the solution to 


Equations (78) and (79) by the requirement that pairs must be 
destroyed at the same rate as they are created in steady state 
(cf. Equation (18) in LZ87). The temperature of this popula¬ 
tion can be determined self-consistently by not allowing any 
net energy transfer between particles and photons via thermal 
Compton scattering (cf. Equations (27) and (28) in LZ87). 

In our implementation, the photon injection term is written 
as 

li-O = IT-0,ej + /ejTlQ.CMB, (80) 


where 

no,ej{x,t) 


GtT TTlgC^ X^ 

Rn{t) gxD ( _ 1 


(81) 


represents the injection of thermal photons from the ejecta 
material with effective temperature Tgff = Teff_com/(C7ej)^^"^ 
(cf. Section 4.3.2). The thermal spectrum fio.CMB of the CMB 
is defined in the same way with an effective temperature of 
Teff.CMB = 2.725 K and it is ‘switched on’ by the function 
/ej once the ejecta material becomes optically thin (cf. Equa¬ 
tion (102)). 

In contrast to LZ87, we additionally include effects of syn¬ 
chrotron radiation (as do Coppi 1992 and Belmont et al. 
2008). Here, we briefly outline the way we include those ef¬ 
fects via the terms Ac.syn and hsyn- In deriving Equation (79), 
we assumed that the particle distribution N{j) can be de¬ 
scribed by the continuity equation 

dN d 

—+ -(7iV)=P(7) + Q(7), (82) 


with dN/dt = 0 in steady state and thus Ac.syn = 
—d{'^N)/dj. Equation (82) is an accurate description in the 
case of synchrotron emission and Compton scattering exclud¬ 
ing the Klein-Nishina limit (as assumed here for simplicity; 
Blumenthal & Gould 1970). Assuming steady state. Equa¬ 
tion (82) can be integrated to give 

^'Tmax 

^(7) = -7c';syn(7) / [Ph') + Qh')] dy, (83) 

where 

TC,syn ( 7 ) =7c(7)+7syn(7) (84) 

is the combined particle cooling rate due to Compton scatter¬ 
ing (excluding the Klein-Nishina limit) and synchrotron emis¬ 
sion, with 

/4 \ r3/47 

icil) = -o-tc (^-7 “ 1 j y n{x)xdx (85) 

(cf. Equation (7) in LZ87) and 

/o„3 TS 

= n{x/Xn)dx ( 86 ) 

hrUnC^ 

(cf. Appendix C, Equation (C4)). The synchrotron losses (86) 
result in a photon source term 

O 1 /-Tmax 

h,,yn{x) = — -^- / N{'^)'R{x/Xc)d-i (87) 

(cf. Appendix C, Equation (C5)). We note that for the station- 
arity assumption leading to Equation (83) to hold in our evo¬ 
lution scheme (Equations (11)-(15)), the timescale for equi¬ 
libration of the particle distribution given by the total cooling 
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timescale | 7 / 7 c,syn| has to be much smaller than any other 
evolution timescale for all energies 7 . We therefore monitor 
this quantity to check the validity of the stationarity assump¬ 
tion during the numerical evolution (see Section 5.4). Finally, 
we note that effects of synchrotron self-absorption have not 
been considered so far and that they can be neglected for our 
purposes. In order to check the validity of this assumption, 
we monitor the optical depth to synchrotron self-absorption, 
which can be approximated by 


^l'syn(^) — 


d 


STTTOgC® 

/ ^( 7 ) 


97 


TZ{x/xc) + f{'r)TZ{x/xc) 


( 88 ) 


d7 


(cf. Appendix C, Equation (C7)). Typically, the nebula is opti¬ 
cally thin to synchrotron self-absorption at energies of interest 
(see Section 5.5). 

Solving the coupled set of Equations (78) and (79) numer¬ 
ically at time t yields a self-consistent set of spectra {n{x), 
N^j), hA{x), np''’(a;), fi^ix), nsyn(x), At^'^(x), Atj^(x), 

Att, hesc{x), P( 7 ), Ac,syn(7)}- It IS important to point out 
that with the modifications described above, these equations 
still conserve energy (see Appendix D). 

The internal energy i?nth of the nebula in terms of the solu¬ 
tion to Equations (78) and (79) can be written as 


Enth{t) = - nRn{t)^meC^ 


n{x, t)x dx 


/ Tniax 

At( 7 ) ( 7 - l)d 7 


(89) 


which is typically dominated by the photon contribution (first 
term on the right-hand side). The pressure of the nebula is 
then given by (cf. also Equation (53)) 


Pnit) = 


^nth(f) ^ 

47ri?3(f) Stt 


(90) 


Moreover, the luminosity of escaping radiation from the neb¬ 
ula is given by 


of kinetic energy deposited into the shocked material by the 
shock. 



Ackin^^^ df 
at 


(94) 


(cf. Equations (65) and (67)). Eor typical sub-relativistic or at 
most mildly relativistic speeds Uej,in of the outer ejecta front 
in Phase I and II (cf. Equation (20)), the initial speed of the 
ejecta shell in Phase III is also at most mildly relativistic. 

However, further acceleration according to (cf. Ap¬ 
pendix B; Equation (11)) 


(95) 


can be significant if the nebula pressure (Equation (90)) is 
high. Here, 7 ej = (1 — denotes the Lorentz factor 

of the ejecta shell and ctej = 47ri?^pn/A7ej is the acceleration 
of the ejecta shell in the frame comoving with the shell. In 
order to define a comoving frame we employ the fact that such 
a freely expanding, spherically symmetric shell at a time t can 
be described as a thin ring of finite extent in a Milne universe. 
The transformation between comoving and lab frame is thus 
achieved via the Milne universe metric. Eor further details, 
we refer to Appendix B. 


Thermal emission — As in Phase I and II, the ejecta matter 
consists of a mixture of photons and ideal gas with adiabatic 
index Tq =4/3, the temperature of which in the comoving 
frame, T^om, can be found by solving 


4 com + 1 Pej^com — a. 


com ' -n 1 ^ 

i ej — 1 m^a 


(96) 


Henceforth qcom refers to the value of a quantity q as measured 
in the frame comoving with the ejecta shell. In Equation (96) 
we have used the fact that i?th/14j = 7?th,com/14j,com (cf. Ap¬ 
pendix B). Furthermore, pq = Mq/Vq^com = Mq/(Vq is the 
density of the ejecta shell, where Vq = ( 4 / 3 ) 7 r(i?gj — R^) is 
its volume in the lab frame and 


3 7|/3ej - arctanh j3q 

~ 2 llPl 


LpwN(a:,f) = -TTRa{t)^meC^nescix,t)x, (91) 

and thus we arrive at the desired expression for the PWN 
source term in Equation (14): 


dE, 


PWN 


dt 


= E 


PWN 


it) = 


LpwNix,t)dx. (92) 


4.3.2. Ejecta properties and emergent radiation 

Shell kinematics — The initial speed 'Uej(fshock,out) of the ejecta 
shell in Phase III is inferred from its kinetic energy i?kin at 

^shock.out- 


j (^shock,out) — C 


1 - 1 


A^kin ( f shock,out) 


MqC^ 


-2 


(93) 


with j3q = Vqjc (cf. Appendix B). The associated effective 
temperature in the comoving frame is given by 




ef¥,com 


{t)^ 


16 


3 ATej,com(t) + 1’ 


where 


ATc) 


1(0 — ^Pej (f) Agj (f)7ej (f) j 


(98) 


(99) 


which gives rise to thermal emission from the outer and inner 
surface of the ejecta shell with luminosities 




and 


Arad,in (i) = 47rii^(f) 


1 


C(07ej(i) 


ot: 


eff.c 




( 100 ) 


( 101 ) 


The total kinetic energy £;kin(f shock,out ) = £^kin,w+i^sh,tot of the 
ejecta matter is the sum of the kinetic energy of the original 
wind material i?kin,w (cf Equation (44)) and the total amount 


respectively (cf. Appendix B; Equation (14)). In Equa¬ 
tion ( 100 ), we have again included a factor /beam = 1 — ^ej/c^ 
to account for relativistic beaming (cf. also Section 4.3.1). 



















Electromagnetic emission from long-lived BNS merger remnants I 


15 


This concludes the discussion of source terms in the main evo¬ 
lution equations (1)-(15) of our model. 


4.3.3. Transition to the optically thin regime 

Once Arej^com (cf. Equation (99)) approaches unity, the 
ejecta material becomes transparent to radiation from the neb¬ 
ula. This non-thermal radiation is not absorbed by the ejecta 
material anymore and becomes directly observable. In order 
to ensure a smooth transition between the optically thick and 
thin regimes, we employ an auxiliary function 


/ejW = 


0 if ATgj (t) > AXgj thres 

^ (f) < ATej,thres 


( 102 ) 


to gradually ‘switch on or off’ terms in the evolution equa¬ 
tions during this transition. Here, 


A -A /l-'t'ej/c 

Arej - Arej,,o™W ^ ^ 


(103) 


approximates the optical depth of the ejecta material as seen 
from the lab frame (cf. Abramowicz et al. 1991), with Arej_com 
defined in Equation (99). We note that this transition function 
and its parameters b and Arej thres are somewhat arbitrary and 
chosen in such a way that they do not influence the numeri¬ 
cal evolution significantly other than guaranteeing a smooth 
transition from the optically thick to the optically thin regime. 
Typically, ATej,thres — few and 6 > 1. 

In particular, as the ejecta material does not absorb neb¬ 
ula radiation anymore, the source term diJpwN/df (cf. Equa¬ 
tion (92)) needs to be removed from the corresponding evo¬ 
lution equation (cf. Equation (14)). Instead, the non-thermal 
emission from the nebula 

imd,nth it) = /ej (t) (104) 

and its associated spectrum given by Equation (91) become 
directly observable. Consequently, the ejecta matter also be¬ 
comes optically thin to photons form the CMB, which are now 
able to diffuse into the nebula (cf. Equations (77) and (80)). 


4.4. Collapse to a black hole 

Our model applies to hypermassive, supramassive, and sta¬ 
ble remnant NSs. If the NS is hypermassive at birth, the ex¬ 
pected lifetime is of the order of fdi ^ 1 s, unless it migrates to 
the supramassive regime through substantial mass ejection on 
shorter timescales. If the NS is supramassive, it can survive 
on much longer spin-down timescales, although (magneto- 
jhydrodynamic instabilities can cause an earlier collapse to 
a black hole. In our model, the time of collapse t = tcoii to a 
black hole is an input parameter that can be adjusted to cover 
all possible scenarios. 

If the NS collapses to a black hole during Phase I, we keep 
evolving Equations (1) and (2), setting Lem dLth,Ns/df to 
zero. As the collapse proceeds on the dynamical timescale, 
which is of the order of milliseconds, we consider the NS 
collapse at f ~ fdr ^ 1 s to be instantaneous as far as the 
numerical evolution is concerned. The time of collapse is 
parametrized by /coii.pi in units of tdr, fcoii = /coii.pifdr- We 
also note that the wind model (cf. Section 4.1.1) has to be ad¬ 
justed by setting M{t') to zero for t' > Loll in Equation (27). 
The resulting EM emission from the system will be predom¬ 
inantly thermal and it will reflect the gradual depletion of the 


energy reservoir of the ejecta material acquired up to the time 
of collapse (see also Paper II). 

In what follows, we consider a collapse occurring in Phase 
III.‘^ In the case that the observed prompt 7 -ray emission of a 
SGRB is associated to this collapse as in the recently proposed 
time-reversal scenario (Ciolfi & Siegel 2015b, a), the resulting 
lightcurves and spectra of our model after t = fcoii correspond 
to the observed afterglow radiation of the SGRB. 

We parametrize fcoii in terms of the spin-down timescale 
icon = /coii^sd- Scenarios for a wide range of values for /coii 
are explored in the companion paper (Paper II). The collapse 
of the NS proceeds on the dynamical timescale, which is of 
the order of milliseconds and it can thus be considered in¬ 
stantaneous at times Loll > Lui.in ^ Is- At times t < 
few X Ld, the typical cooling timescales |7/7c,syn| of the non- 
thermal particles in the nebula (cf. Equation (84)) are orders of 
magnitude smaller than any other evolution timescale of our 
model thanks to the relatively strong magnetic field (cf. Equa¬ 
tion ( 86 ); see also Section 5.4 and Section 4.2 of Paper II). 
Eurthermore, due to this efficient cooling their total number 
density Ajot = / ^( 7 , t) d 7 Ajh.tot is orders of magnitude 
smaller than the total number density Ajh,tot of thermalized 
particles in the nebula. We can therefore assume that the non- 
thermal particle population is instantaneously thermalized at 
t = fcoii- Eurthermore, noting that the optical depth to Thom¬ 
son scattering Atj is given by Att = (Txi?nAth,tot, we can 
conclude that Att is not affected by this instantaneous ther- 
malization. Pair annihilation becomes increasingly unlikely 
as the nebula further expands, and can partially be compen¬ 
sated by pair creation induced by photons in the high energy 
tail of the photon spectrum above the pair creation threshold. 
As a consequence, the total number of particles is roughly 
conserved, and the evolution of the Thomson scattering depth 
for t > Loll is approximately given by 


ATj{t) 


ArT(Loii) 


-Rn(Loll) 

Rnit) 


(105) 


Due to efficient particle cooling prior to f = fcoii (see above), 
the energy budget of those photons above the pair creation 
threshold is typically orders of magnitude smaller than the 
combined photon energy below the pair creation threshold at 
fcoii- Furthermore, we expect photons to interact with par¬ 
ticles mostly via Thomson scattering, as photons that could 
Compton-scatter off thermal particles are small in number. 
Noting that Thomson scattering is elastic, we can deduce that 
for our purposes, the spectral shape of the non-thermal pho¬ 
ton spectrum att = Loll remains frozen thereafter as photons 
diffuse out of the nebula on the diffusion timescale 

Liff,n(i) “ (106) 

Here, we have again added the light crossing time of the neb¬ 
ula to account for a transition to the optically thin regime. 
The luminosity of the photons diffusing out of the nebula is 
approximately given by 


dLpwN 

df 


— .^PWN (f) 


(107) 


= {!-[!- /ej(t)]^(f)}/beam(f) " , y 

^diff,n \^) 

Here we do not discuss the possibility of a collapse during shock prop¬ 
agation (Phase II), as this phase is typically very short compared to Phase I 
and III, and a collapse is thus unlikely to occur. 
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where, in analogy to Equation (57) the prefactor takes into 
account the combined effects of relativistic beaming and a 
non-zero albedo of the surrounding ejecta material.^ The total 
energy budget of the nebula thus evolves according to 


wind density profiles at later times. However, as the shock 
sweeps up all the material into a thin shell in Phase II, the 
exact distribution of matter is not important and does not in¬ 
fluence the lightcurves and spectra of our model significantly. 


,, — ~-^^PWN(i) + Aad,in(f)j (108) 

at 

where Lrad.in is the thermal radiation from the surrounding 
ejecta shell (see Equation (101)), and we can write the non- 
thermal radiation leaving the system as 

Lrad.athiXjt) = /ej (f)EpWN(a;, fcoll) y——• (109) 

-EpwN(tcoll) 

In conclusion, after a collapse to a black hole, we evolve the 
system of evolution equations in Phase III assuming conser¬ 
vation of magnetic energy (dEs/dt = 0), employing Equa¬ 
tion (107) in Equation (14), and extending the set of Equa¬ 
tions (11)-(15) by Equation (108) to evolve the nebula prop¬ 
erties instead of solving Equations (78) and (79). 

5. NUMERICAL PROCEDURE 

In this section, we discuss several aspects related to the 
numerical evolution of the model presented in the previous 
sections. After briefly outlining the overall numerical pro¬ 
cedure to integrate the evolution equations (1)-(15) in Sec¬ 
tions 5.1, 5.2, and 5.3, we identify and define important 
timescales to be monitored during the numerical evolution 
(Section 5.4), discuss how to monitor the validity of neglect¬ 
ing synchrotron self-absorption (Section 5.5), discuss how to 
overcome problems related to stiffness in Equations (1)-(15) 
(Section 5.6), and Anally discuss how to numerically compute 
the lightcurves as seen by a distant observer including rela¬ 
tivistic effects such as relativistic beaming, the time-of-flight 
effect, and the relativistic Doppler effect (Section 5.7). 

5.1. Hydrodynamic wind evolution 

Prior to evolving the main evolution equations (1)-(15), the 
evolution of the background fluid (the baryonic wind emit¬ 
ted during Phase I; see Section 4.1.1) has to be computed 
numerically. Eirst, at any given time t the function t{r, t) is 
constructed solving Equation (28) numerically for every ra¬ 
dius r. Knowing i{r, t) immediately yields the velocity pro¬ 
file Vw{r,t) via Equation (30) and the mass profile mv,{r,t) 
via Equation (27). By finite differencing the mass func¬ 
tion mv,{r,t) the corresponding density profile /Ow,t(i’) is ob¬ 
tained (cf. Equation (29)). We typically employ a spatial grid 
ranging between r = R^in — 30 km and i?ej(f) (cf Equa¬ 
tion ( 20 )), which is logarithmically spaced with a fixed num¬ 
ber of points per decade. Automatic mesh refinement (where 
needed) is used to resolve the profiles, which become increas¬ 
ingly ‘sharp’ as the material moves outward to larger length 
scales. By employing this wind model the hydrodynamic evo¬ 
lution of the wind material has essentially been reduced to a 
numerical root finding problem. 

The free parameters of the model, Min, fdr, o'M, and Wejjn 
are listed in Table 2. We note that numerical relativity simu¬ 
lations can be employed to directly determine Min and Uej.in- 
Eurthermore, they can provide estimates on fdi. The only un¬ 
constrained parameter is ctm, which controls the shape of the 

^ Prior to collapse in Phase III, the quasi-stationarity assumption in the 
present implementation of the radiative processes of the PWN (Section 4.3.1) 
prevent us from introducing these effects (see also Section 6). 


5.2. PWN balance equations 

At any time t in Phase III, the balance equations (78) and 
(79) for the photon and non-thermal particle distribution in 
the PWN have to be solved. These distributions then spec¬ 
ify, e.g., the source term di?pwN/df (cf Equation (92)) re¬ 
quired in Equation (14). We solve the coupled set of highly 
non-linear, integro-differential equations (78)-(79) in analogy 
to the multi-step, iterative method outlined in Appendix B of 
LZ87. Only Step 1 of their method differs slightly from our 
Step 1 due to the modifications we introduced to the physi¬ 
cal description of the pair plasma (cf. Section 4.3.1). There¬ 
fore, we rediscuss this step here for completeness and refer 
the reader to LZ87 for the other Steps. 

Step 1 in solving Equations (78)-(79) consists of neglect¬ 
ing thermal Comptonization (i.e., setting nj = 0 in Equa¬ 
tion (78)) and proceeds as follows. 


(i) Given an approximate photon spectrum n{x) com¬ 
pute the corresponding particle distribution N ( 7 ) using 
Equation (83), where the integral over P{'y) is com¬ 
puted using Equation (13) of LZ87. 


(ii) Compute 

R 

A{x) = -^{nQ+riA + + flsyn) 

using W( 7 ) from (i). Equation (80), Equa¬ 
tions (7),(8),(11),(15)-(18) of Svensson (1983), 
Equation (9) of LZ87, and Equation (87), 


B{x) 


1 

1 + lTKN{x)f{x) 




using Equations (21) and (10) of LZ87 together with 
N{'y) from (i). Eurthermore, compute 


C{x) = 0.2i?nCrT- 

X 

as well as D{x) = A[l/x), E{x) = Biljx), and 
F{x) = C{l/x). We note that, with these defi¬ 
nitions, all coefficient functions A,B,C,D,E,F are 
non-negative. 


(hi) Writing Equation (78) as n{x) = A{x)/[B{x) + 
C(a;)n(l/x)] and n{l/x) = D{x)/[E{x) + F{x)n{x)] 
leads to a quadratic equation for n(x). 


BFn^{x) + {BE + CD- AF)n{x) - AE = 0, 


which has only one physical root, 
{BE + CD- AF) 


i{x) = - 


-f 


2BF 

y/{BE + CD- AFY E ABFAE 
2BF ■ 


( 110 ) 


Equation (110) is now used to update n{x). 

Initializing n{x) by (i?n/c)[ho,ej + /eji^o.CMfi], we iterate 
on Steps (i)-(iii) until for successive iterations j — 1 and j we 
have 

e = ||?T.(x)j —1 Umax ^ 


( 111 ) 
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Here ||u(a;)||max = for a function u{x) on the 

dimensionless photon energy domain X = [tCmin, a^max] and 
typically Ctoi < 0.01. As 7niax (cf. Equation (76)) is the max¬ 
imum attainable photon energy, we set ccmax = 7max- In order 
to cover the frequency range down to the radio band, we typi¬ 
cally set Xmin = 10“^®. Both the particle energy and the pho¬ 
ton energy domains are logarithmically spaced with typically 
15 points per decade. 

5.3. Evolving the main model 

After having computed the hydrodynamic evolution of the 
baryonic wind (Section 5.1), its density profile Pw,t(T) at any 
point in time and space is known and the main equations of 
our model (1)-(15) can be integrated. The time grid T = 
[fmin, fmax] for this evolution is logarithmically spaced with 
typically 50 points per decade. Additional mesh refinement is 
used in order to properly resolve the propagation of the rela¬ 
tivistic shock front through the ejecta material in Phase II. In 
our model, we associate the time of merger of the BNS system 
with t — 0, but start the numerical evolution from an appropri¬ 
ate final checkpoint of a numerical relativity simulation, cor¬ 
responding to tens of milliseconds after merger, and evolve 
it over the timescales of interest until, e.g., f^ax lO^s. It 
is important to point out that numerical relativity simulations 
of BNS mergers can be employed to determine most of the 
input parameters of our model, which are listed in Table 2. 
Given those initial parameters read off or estimated from a 
simulation, the following evolution according to our model is 
a self-consistent prediction based on these initial conditions. 

5.3.1. Phase I 

In order to evolve Equations ( 1 ) and (2) from time tlot+Xt 
we need to compute the source terms on the right-hand sides. 
This can be done in the following way: 

(i) Compute Vv,{Rej{t),t) using the precomputed baryonic 
wind model (Section 5.1) or Equation (20). 

(ii) Employ the previously generated wind profiles 
TOw(t, f) and /Ow(f) to compute the optical depth of 
the ejecta material (Equation (36)) and to numerically 
solve for the temperature of the ejecta material using 
Equation (39). This yields the source term Erad(f) 
(cf. Equations (34) and (35)). 

(iii) Compute the remaining two source terms using Equa¬ 
tions (32) and (33). 

This phase ends when t = fpui.in (Equation (40)) and Phase 
II starts. If the NS collapses to a black hole during Phase 
I, we keep evolving Equations (1) and (2), setting Lem and 
dLth,Ns/df to zero (cf Section 4.4). 

5.3.2. Phase II 

At f = fpui,in, the initial rotational energy of the pulsar is 
computed using Equation (42) and its initial spin period is 
then inferred from Equation (46). This allows us to com¬ 
pute the initial spin-down luminosity (Equation (47)) and the 
spin-down timescale (Equation (49)), which together define 
the spin-down luminosity at later times (Equation (48)). 

In order to evolve the evolution equations (3)-(10) from 
time t to t + At, one can proceed as follows: 

(i) Compute using the precomputed baryonic 

wind model (Section 5.1) or Equation (20). 


Table 2 

Model input parameters. Most of these parameters can be extracted from (or 
at least estimated/constrained using) numerical relativity simulations of 
BNS mergers. 


Parameter Description 

Min initial mass-loss rate of the NS (cf. Section 4.1.1) 
fjr timescale for removal of differential rotation from the NS 
(cf. Section 4.1.1) 

<JM of fdr to the timescale for decrease of the mass-loss rate 

(cf. Section 4.1.1) 

r;ej,in initial expansion speed of the baryonic ejecta material 
(cf. Section 4.1.1) 

B magnetic held strength in the outer layers of the NS (cf. Equa¬ 

tion (32)) 

r)B dipolar magnetic held strength of the pulsar in units of B 

(cf. Equation (41)) 

^rot,NS,in initial rotational energy of the NS (cf. Equation (42)) 

Pc initial central spin period of the NS (cf. Equation (32)) 

Re equatorial radius of the NS (cf. Equation (32)) 

MNs,in initial mass of the NS (cf. Equation (45)) 

7pui moment of inertia of the pulsar (cf. Equation (46)) 

Cej.NS.in initial specihc internal energy of the material ejected from the 
NS surface (cf. Equation (33)) 

K, opacity of the ejecta material (cf. Equations (36), (71), (72), 
and (99)) 

A{u., t) frequency and time-dependent albedo of the ejecta shell in 
Phase II and III (cf. Sections 4.2.2 and 4.4) 
tiy neutrino-cooling timescale (cf. Equation (33)) 

fraction of the total pulsar wind power injected as magnetic 
energy per unit time into the PWN (cf. Equations (10), (15), 
and Sections 4.2.2, 4.3.1). 

?77S efficiency of converting pulsar wind power into random ki¬ 
netic energy of accelerated particles in the PWN (cf. Equa¬ 
tions (9), (76) and Sections 4.2.2, 4.3.1) 

7max maximum Lorentz factor for non-thermal particle injection 
into the PWN (cf. Section 4.3.1, Equation (76)) 

Ee power-law index of the non-thermal spectrum for particle in¬ 
jection into the PWN (cf. Section 4.3.1, Equation (76)) 

/coll (only in the collapse scenario. Section 4.4) parameter spec¬ 
ifying the time of collapse of the NS in units of the spin- 
down timescale (collapse during Phase III) or in units of 
(“/coii.Pi”’ collapse during Phase I) 


(ii) Compute dLth.voi and fsh(i) from Equations (62) and 
(58), respectively. Use v^hit) and the wind profiles to 
evaluate shock heating according to Equation (63). 

(iii) Compute the source terms Lrad.pub and 

dLpwN/df using Equations (48), (51), and (57). 

(iv) Eind the temperatures of the shocked and unshocked 
ejecta parts by using the precomputed wind profiles 
mw{r, f), Pw(f) and solving Equation (69). This yields 
the optical depths (Equations (71), (72)) and luminosi¬ 
ties (73) and (74). 

(v) Evolve all quantities to f -I- At, except for R„ and 

(vi) Dehne 

Enlh{t + At) = Enthit) H- —At (112) 

Rr,it + At) = + At) (113) 

Enth(t + Af) 

and iterate until both quantities have converged to some 
desired accuracy. The second condition ensures pres¬ 
sure equilibrium between the nebula and the shocked 
ejecta layer (cf. Equation (61)). 

Phase II ends when i?sh(f -I- At) > i?ej(f + At). 
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5.3.3. Phase III 

Once the shock has reached the outer ejecta layers. Phase 
III begins. The initial speed of the shocked ejecta layer is 
calculated from Equation (93). Evolving Equations (11)-(15) 
from time t to t + At can be accomplished by the following 
steps: 

(i) Compute L^d and Erad.pui using Equations (48) and (51). 

(ii) Compute the acceleration of the ejecta shell according 
to Equation (95). 

(iii) Eind the temperature of the ejecta shell by solving 
Equation (96). This yields the optical depth (Equa¬ 
tions (99)) and the luminosities Lrad and Trad,in (Equa¬ 
tions ( 100 ) and ( 101 )). 

(iv) Solve Equations (78) and (79) as described in Sec¬ 
tion 5.2 to hnd the source term d£^pwN/df. 

In the case the NS collapses to a black hole (see Sec¬ 
tion 4.4), Step (i) is omitted and Equation (108) is added as 
another evolution equation to Equations (11)-(15). Eurther- 
more. Step (iv) is replaced by 

(iv)’ Eind the source term d£'pwN/df using Equation (107). 

This concludes our discussion of the overall procedure to 
evolve Equations (1)-(15). 


5.4. Timescales 


As our evolution model in Phase III is built upon the as¬ 
sumption of quasi-stationarity as far as radiative processes 
in the PWN are concerned, we need to monitor several 
timescales during the numerical evolution in order to assess 
the validity of this assumption. Eor these diagnostic purposes, 
we dehne the following timescales in Phase III: 


A(f) = p 

^ph 

Tc{l,t) = 


7C,syn 

Rn 

C 


(114) 

(115) 

(116) 
(117) 


where 7 denotes the Lorentz factor of a particle in the PWN 
and Ig, /ph, and 7 c,syn are dehned by Equations (76), (77), and 
(84), respectively. 

In order for the stationarity assumption regarding the parti¬ 
cle distribution to hold, which allowed us to integrate Equa¬ 
tion (82) to obtain Equation (83), the timescale for equilibra¬ 
tion of the particle distribution given by the cooling timescale 
Tc has to be much smaller than any other timescale involved 
in the problem. In particular, it has to be smaller than the 
timescale for change of the photon distribution. Thus we have 
the requirements that 


Tc{l,t)<.Tpd{t) Vf, 7 , (118) 

Vt, 7 . (119) 


Here, we have assumed that the timescales for change of the 
photon and particle distributions inside the nebula are approx¬ 
imately given by the timescales Tph and Te for change of the 


injected spectra. As we shall discuss in the companion pa¬ 
per (Paper II), we typically hnd that these conditions are very 
well satished across the entire parameter space in absence of 
signihcant acceleration of the ejecta shell after t = fshock,out 
and except for very late times t > 10 ^ s when the nebula has 
grown in size, the magnetic held strength has decreased, and 
synchrotron cooling has become less efficient. However, we 
are not interested in the evolution at these late times. 

Eor the stationarity assumption regarding the photon spec¬ 
trum to hold (i.e., for Equation (78) to hold), the timescale for 
the photon spectrum to change has to be much larger than the 
timescale for equilibration inside the nebula. As the positronic 
plasma inside the PWN is relativistic, its sound speed is close 
to the speed of light. We can thus assume that the nebula 
plasma adjusts to changing exterior conditions essentially on 
the timescale t\. Therefore, one also has the requirement that 

'ri(t) < Tph(f) yt. ( 120 ) 

As long as further acceleration of the ejecta shell after t = 
fshock.out is not signihcant, we typically hnd that this condition 
is well satished during the numerical evolution as we shall 
discuss in more detail in the companion paper (Paper II). 

5.5. Synchrotron self-absorption 

In our model of the PWN in Phase III, we have neglected 
synchrotron self-absorption. In order to check the validity of 
this assumption during the numerical evolution of our model, 
we monitor the dimensionless photon energy a:ATsy„=i(i) de¬ 
hned by 

Arsyn(a:AAy„=i(t),t) = 1, (121) 

where AT^ya{x,t) is the approximate optical depth of the 
nebula to synchrotron self-absorption at time t (cf. Equa¬ 
tion ( 88 )). This dimensionless energy separates the optically 
thick and thin regimes of the photon spectrum and lies well 
below X-ray frequencies (typically by many orders of magni¬ 
tude) as we shall discuss in the companion paper (Paper II). 
Therefore, the part of the spectrum that we are mostly inter¬ 
ested in, i.e., at X-ray and 7 -ray energies, is unaffected by 
effects of synchrotron self-absorption. 

5.6. Stiffness problem 

At some point in the evolution (in Phase III), the photon 
diffusion timescale of the ejecta shell 

~ (Arej + 1) (122) 

c 

becomes comparable to the temporal resolution At and Equa¬ 
tion (14) becomes stiff Here, At^ denotes the optical depth 
as seen from the lab frame (cf. Equation (103)). This problem 
can be noticed in the following way. Assuming ATej ^om 3> 
1 , a radiation dominated gas, such that T^om — E±jaVej 
(cf. Equation (96)), and using Equations (98) and (100) it is 
straightforward to show that the thermal emission from the 
ejecta shell scales as L^-^d oc £’th/idiff,ej- The photon diffusion 
timescale itself scales as fdiff.ej oc (cf Equation (99)). 
Hence, as i?ej is monotonically increasing with time, at some 
point fdiff.ej ~ At and Equation (14) becomes stiff. 

Eortunately, there is an easy way to bypass this problem. 
As Lrad oc i?th/fdiff.ej, this term will eventually deplete the 
energy reservoir E± of the ejecta shell until the shell enters 
an asymptotic regime in which the amount of injected energy 
equals the amount of emitted energy on the photon diffusion 
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Figure 3. Evolution of the internal energy -Eth of the ejecta shell in Phase 
III according to a run using the evolution scheme described in this section 
(Run 1) and a run using additional mesh refinement to evolve Equation (14) 
up to late times (Run 2). The dashed red line indicates the time of transition 
t = tres to the modified evolution scheme for Run 1. 


timescale. In other words, this regime is defined by setting 
d£^th/d^ = 0 in Equation (14) on a timescale t ~ idiff.ej- At 
any time in this asymptotic regime, the total internal energy 
Eih of the ejecta shell is then given by 

i?th(i) = [1 - /ej(t)]^^Uff,ej(f) (123) 

and a separate evolution equation for E(h is not needed any¬ 
more. 

In practice, the grid spacing At has to be adjusted such that 
this asymptotic regime is reached well before At « fdiff.ej- 
Once this is achieved, we typically switch to the asymp¬ 
totic regime during the numerical evolution at ires, defined 
by fdiff,ej(ires) = 2 At. By definition of this regime, and 
Arad.in scale as d£^pwN/df. Therefore, assuming that the ther¬ 
mal photons from the ejecta shell dominate the photon injec¬ 
tion Zph of the PWN (cf. Equation (77)), Zph is expected to 
scale as 1^. As the PWN conserves energy (total luminosity of 
injection equals the luminosity of emitted radiation; see Ap¬ 
pendix D) we hence conclude that 


d£^pwN 

dt 


tN (Asd L rad,pul) ■ 


(124) 


In order to avoid numerical runaway instabilities due to the 
fact that di?pwN/df in Equation (123) depends on the solu¬ 
tion to Equations (78) and (79), which, in turn, depends on 
the injected photons, i.e., on di?pwN/d<, we employ Equa¬ 
tion (124) in Equation (123) and calibrate the prefactor at ires 
to ensure continuity over t = fres- This calibration, however, 
is only valid for t > t^es if also the dynamics of the shell have 
reached an asymptotic regime, i.e., if no further acceleration 
occurs such that, e.g., the beaming factor does not depend on 
time. In practice, a grid spacing can always be set up to satisfy 
this constraint as well. 

Figure 3 compares the internal energy i?th of a typical 
model case evolved using the scheme outlined above (Run 
1) with the corresponding result for a simulation in which 
a power-law refinement At oc fdiff.ej oc was imple¬ 

mented (Run 2). For the latter run, the refinement guaranteed 


At <C fdiff.ej at all times and thus allowed us to integrate Equa¬ 
tion (14) for all times. The agreement between the two runs 
is remarkably good and indicates that the scheme described 
above well captures the evolution of the system at later times 
than fjes- Remaining discrepancies at late times t ^ 10^ s are 
due to numerical errors as the ever decreasing time step in 
Run 2 results in increasing accuracy for the time integrations 
of Equations (11)-(15). Finally, we note that runs with power- 
law refinement are computationally very expensive, such that 
comparisons like the one presented here cannot be carried out 
routinely and we will thus not add this comparison to the set 
of routine checks (see Sections 5.4 and 5.5) that we shall dis¬ 
cuss in more detail in the companion paper (Paper II). 


5.7. Observer lightcurve reconstruction 

This section is devoted to a discussion of how the lightcurve 
and spectra as seen by a distant observer can be reconstructed, 
including relativistic effects such as the relativistic Doppler 
effect, the time-of-flight effect, and relativistic beaming. As 
the effective temperature of blackbody radiation is altered by 
the relativistic Doppler factor and enters the luminosity to the 
fourth power, even mildly relativistic ejecta shell velocities of 
Uej ~ 0.1c can have a significant influence on the lightcurve. 
Furthermore, as the PWN and the ejecta shell expand to large 
radii, radiation reaching the observer from different locations 
on the surface of the expanding sphere can have appreciable 
delays. These delays are particularly important in dynamical 
situations, in which the ejecta shell is accelerated or in which 
the NS collapses to a black hole and induces an abrupt change 
in the radiation escaping from the system. Finally, relativis¬ 
tic beaming is influential in the sense that it normalizes the 
total luminosity by selecting an effective surface area of the 
expanding sphere that emits in the direction of the observer. 
In so doing, it also affects the maximum delay concerning the 
time-of-flight effect and must therefore be taken into account. 
Consequently, it is important to consider these relativistic ef¬ 
fects when predicting observer lightcurves and spectra with 
our model. 

We numerically compute the lightcurves and spectra as 
seen by a distant observer including the relativistic effects de¬ 
scribed above by decomposing the radiation originating from 
the surface area of the half-sphere facing the observer into in¬ 
dividual energy packages released at time t during a time At 
from ring-shaped portions of surface area of equal distance to 
the observer. These packages are emitted at time t and then 
recollected by the observer according to their time of flight 
and their Doppler-shifted frequency. To this end, we define 
a coordinate system centered about the NS, with its z-axis in 
direction of the remote observer and 9 denoting the polar an¬ 
gle. At any time t let us consider two half spheres facing the 
observer, defined hy 9 G [0, 7r/2] and radii i?ej(f) and Rn{t), 
respectively. Due to relativistic beaming, only radiation orig¬ 
inating from the regions 9 G [0,0rnax(i)] will eventually reach 
the observer (Rees 1966), where 


0max(i) = arccos 



(125) 


This area is then further decomposed into rings of constant 
9, defined by a grid 0 = [0,..., 0^,..., 0max(i)]. with ng 
points and spacing A9k that is equidistant in the limb angle 
fi = cos 9. 

Moreover, given a time grid T = [fmin,..., ..., fmax] for 

the evolution of Equations (1)-(15) appropriate to properly 
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resolve the shock dynamics during Phase II and to overcome 
the stiffness problem (cf. Section 5.6), we define a time grid 
for the observer T' = 
by 


f' 

If • • • 5 5 


! ^max ~ fmax] 


f' — f' 

''0 ''mm) 


f' - f' 
+ l “ S' 


cne 


(126) 

(127) 


where Nf < 1. In order for the observer to resolve the time- 
of-flight effect, i.e., to distinguish between the various energy 
contributions from different locations on the emitting surface 
area, the grid spacing Af' at time t' has to be smaller than the 
light crossing time of the half sphere divided by the number 
of individual rings. Since 


AS = 

cne cne 


(128) 


the grid spacing defined by Equation (127) guarantees that the 
time-of-flight effect is properly resolved. 

We also define a receiver energy grid E'{T',X') for the 
remote observer, where X' denotes the grid in dimensionless 
energy x' = hn 'which is typically similarly spaced as 
the corresponding one used in the evolution of the main model 
(cf. Section 5.2). The energy contribution of a ring of width 
dO at 0 to the total energy emission during a time dt per unit 
frequency as measured by the distant observer is given by 


= AttI'{ u')da cosOcomdt, (129) 

du' 

where I'{v') is the specific intensity as seen by the observer, 
dcr = 27ri?gj/^ sin OdO is the surface area of the ring, and 6>com 
given by 

/ T , „ sin0 

tan ^com = \/1 - c^ (130) 

V *=1' COS^-Uej/c 

defines the direction of the observer as seen in the local co¬ 
moving frame of the expanding surface. We note that the ef¬ 
fective emitting surface area dcr cos 0com approaches zero as 
0 reaches ^max- We have already multiplied by drr in Equa¬ 
tion (129) to account for an isotropic source. 

According to the relativistic Doppler effect applied to a 
photon of frequency v emitted from the surface at a point 
specified by an angle 0 in direction of the remote observer. 


77.^ = ^ 

V X 



1 — (Wj/c) COS0 


(131) 


Noting that I jv^ is Lorentz invariant, where I is the specific 
intensity, we have I' = D^I. Therefore, for the blackbody 
emission of the ejecta material. 


W) = 




exp 


hv' 


77A:B7’eff,c. 


- 1 


(132) 


where Tss^com is given by Equation (98). We note that the 
effective temperature as seen by the distant observer is thus 
^eff = 77Teff_com- Since the associated luminosity scales as 
L' cx small deviations of D from unity can already 

have an appreciable influence on the observer lightcurve. Fur¬ 
thermore, as long as the nebula is optically thick, i.e., Att > 


1, the escaping photons originate from the outer surface layer 
at r = i?n and we can set 


= D 


3 7/rad,nth(a^ /77, t) 


47ri?2 


(133) 


for the non-thermal radiation, with Trad,nth (a;, t) being the non- 
thermal luminosity of the nebula (cf. Equations (104) and 
(109)). When beaming is already encoded in Lrad,nth(a;, t) it¬ 
self (cf. Section 4.4), the right-hand side of Equation (133) 
has to be divided by /beam to solely take the geometric time- 
of-flight effect into account (through the subdivision into in¬ 
dividual rings up to 6* = ^max)- 
When the nebula is optically thin, i.e., Atj < 1, the escap¬ 
ing photons are emitted uniformly and isotropically through¬ 
out the volume of the nebula. In this case, we slice the nebula 
volume 14 into spherical shells of radius r G R = [iimin; Rn] 
and thickness dr and subdivide each shell into rings of con¬ 
stant 0voi G ©voi = [0) tt] and width d^voi- The contribution of 
such volume elements dV = dcrdr = 27rr^ sin ^yoid^voidr to 
the total energy emission during a time df per unit frequency 
as measured by the distant observer can then be written as 


dE' 

dx' 


D^LradMhix'/D, t) l^dt. 

Tn 


(134) 


Here, D is defined as in Equation (131) in terms of the local 
velocity v{r) = v^^r/Ra- 

With these definitions and expressions at hand, we compute 
the radiation emitted at time t during a time At as received by 
the remote observer according to the following steps: 


(i) First, compute the arrival times of the energy packages 
for all emitting rings and spherical segments, i.e., for 
all 0 G O, 0m\ G ©voL and r G R: 

Cth(0) = —+7-^^cos0 (135) 

’ c c 

= — (136) 

c c 

Cnthvol(6'vol,r’) = + t-- COs6»vol (137) 

’ ’ C C 

This definition of arrival times ensures that a photon 
emitted at f = tmm from the outer surface of the bary- 
onic wind at i?ej(7min) = Rmm and 0 = 0 reaches the 
observer at t' = 


(ii) Second, compute the energy spectrum Ai?4(©)2f') 
that is generated by thermal radiation from the surface 
of the ejecta matter using Equations (129) and (132): 

AE;^i0k,iyl)=87r^Rll'M) 

X sin6»fc cos[6>com(6'fe)]A6>fcAf. 


(iii) In analogy to the thermal radiation in Step (ii), use 
Equations (129), (133), and (134) to compute the en¬ 
ergy contributions generated by the non-thermal radi¬ 
ation from the PWN once the ejecta shell has become 
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optically thin: 

Ai?U,surf(^fc,^;) = (1 - fn)87T^RllM) 

X sin4 cos[6>com(6'fc)]A0feAt, 

3 

AEjjd, Tm) — L/rad,nth{Xi / D) 

j,2 

X —^ sin ^voijfc A^voi,fc Ar772AE 

-Tin 

The function fn{t) is defined in terms of Atj in anal¬ 
ogy to /ej in Equation (102). It is used here to guaran¬ 
tee a smooth transition between the optically thick and 
thin regimes of the nebula, corresponding to surface and 
volume nebula emission, respectively. 

(iv) Use the mapping of angles to arrival times from Step 
(i) to put the energy packages from Steps (ii) and (iii) 
into the correct time and frequency bins of the observer 
energy grid E'{T',X'). 

Once the evolution of Equations (1)-(15) has been accom¬ 
plished up to f = ftnax and all energy packages have been sent 
and received, the energy array of the observer, E'{T', X'), 
may be divided by the time steps Af' to obtain the corre¬ 
sponding observer luminosity 

Lobsit'^,x[) = E'{t'^,x'i)/At'^. (138) 

This quantity can then be used to compute detailed predictions 
for the luminosity as seen by a remote observer in specific 
wavelength bands, such as in the 7 -ray, X-ray, UV, optical, 
and radio bands (see Paper II). 

6. DISCUSSION AND CONCLUSION 

In this paper, we have presented a dynamical model to de¬ 
scribe the post-merger evolution of a BNS system and its 
EM emission. Our model assumes that the merger of two 
NSs leads to the formation of a (hypermassive, supramas- 
sive, or stable) NS which does not collapse to a black hole 
on timescales of at least tens of milliseconds after merger. As 
we have argued (cf. Section 1), such a long-lived NS is a very 
likely possibility, such that the model should be applicable to 
a large fraction of BNS merger events. In contrast to a black 
hole promptly formed after merger and surrounded by a short¬ 
lived accretion disk, such long-lived objects can power long- 
lasting EM emission from 7 -ray to radio energies that might 
be responsible for at least part of the observed long-lasting 
afterglow emission observed in many SGRBs. We refer to 
Paper II for a detailed account on modeling X-ray afterglow 
lightcurves with our model in the context of SGRBs. More¬ 
over, this model also represents an important tool to identify 
EM counterparts associated with the GW signal of the inspi¬ 
ral and merger of a BNS system (see Paper II for a detailed 
discussion). The identification of such EM counterparts is es¬ 
sential for performing joint EM and GW observations. Joint 
observations involving EM counterparts of the kind discussed 
here (long-lasting and highly isotropic) can confirm the as¬ 
sociation of the GW signal with a BNS merger (i.e., distin¬ 
guish form the NS-BH binary merger). With the advanced 
LIGO/Virgo detector network starting its first science run later 
this year, such multimessenger astronomy will turn into excit¬ 
ing reality in the very near future. 

Our model is formulated in terms of a set of highly coupled 
differential equations, which provide a self-consistent evolu¬ 
tion of the post-merger system and its EM emission given 


some initial data. Such initial data can be extracted from a 
numerical relativity simulation of the merger and early post¬ 
merger phase at a few to tens of milliseconds after the merger, 
once a roughly axisymmetric state has been reached. Our 
model allows us to evolve the post-merger system over time 
and lengthscales inaccessible to numerical relativity simula¬ 
tions, typically up to ~ 10^ s after merger. It thus bridges the 
gap between numerical relativity simulations of the merger 
process and the timescales of interest for SGRB afterglow ra¬ 
diation. The model evolves the system through three main 
evolutionary phases: an early baryonic wind phase (Phase 
I), a pulsar wind shock phase (Phase II), and a PWN phase 
(Phase III). Eurthermore, the possibility of collapse to a black 
hole during any of the three phases is accounted for. Our 
model links the evolution of the central engine directly to 
the observed afterglow radiation, taking into account relativis¬ 
tic dynamics as well as an accurate reconstruction of the ob¬ 
server lightcurve including relativistic beaming, the relativis¬ 
tic Doppler and the time-of-flight effect. 

Prompt SGRB emission and time-reversal scenario — Our model 
makes no assumption on how and when the prompt 7 -ray 
emission of the SGRB itself is produced. It can accommodate 
both the standard scenario (SGRB at the time of merger) and 
the recently proposed time-reversal scenario (SGRB at the 
time of collapse of the remnant NS; Ciolh & Siegel 2015b, a). 
It can obviously also accommodate the case in which no rela¬ 
tivistic jet and thus no SGRB is produced at all. 

In the context of magnetar models for SGRBs, the prompt 
emission is assumed to be generated by an accretion pow¬ 
ered relativistic jet emerging from a NS-torus system shortly 
(^ ms) after the time of merger (e.g., Metzger et al. 2008; 
Bucciantini et al. 2012; Gompertz et al. 2014; Metzger & 
Piro 2014; Gao et al. 2015). However, a generic feature of a 
newly-bom NS after a BNS merger will be strong baryon pol¬ 
lution in its vicinity due to dynamical ejecta from the merger 
process and neutrino and magnetically driven winds from its 
surface or possibly from an accretion disk (Hotokezaka et al. 
2013; Oechslin et al. 2007; Bauswein et al. 2013; Kastaun & 
Galeazzi 2015; Dessart et al. 2009; Siegel et al. 2014; Met¬ 
zger & Eernandez 2014; see Section 4.1). Such baryon pollu¬ 
tion can not only choke jets (Nagakura et al. 2014; Murguia- 
Berthier et al. 2014), but it is likely to even prevent the gener¬ 
ation of any relativistic outflow at all. We note that numerical 
simulations of BNS mergers that lead to the formation of a 
remnant NS have not found indications for the generation of 
a relativistic jet (Giacomazzo & Perna 2013). 

In the time-reversal scenario, the SGRB is generated at the 
time of collapse of the supramassive remnant NS. In this case, 
the baryon-free PWN surrounding the collapsing NS does not 
threaten the formation of an accretion powered jet from the 
remaining BH-torus system. Margalit et al. (2015) have re¬ 
cently argued that the formation of an accretion disk following 
the collapse of the supramassive NS is rather unlikely. How¬ 
ever, numerical simulations will be needed to further investi¬ 
gate this issue. 

As there are still many open questions related to the for¬ 
mation of the SGRB prompt emission in BNS mergers, our 
model is designed to be general in the sense that it does not 
make assumptions on the generation of the prompt emission. 
In either scenario, it provides detailed predictions for the in¬ 
trinsic afterglow emission emerging from the post-merger sys¬ 
tem that can be compared to observations (see Paper II). In the 
time-reversal scenario, it additionally predicts the EM emis- 
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sion prior to the SGRB itself (see Paper II). If such emission 
is found, it would represent strong evidence in favor of the 
time-reversal scenario (see Section 7 of Paper II for a more 
detailed discussion and implications for joint GW and EM 
observations). 

Comparison to earlier work — Other authors, e.g., Yu et al. 
(2013) and Metzger & Piro (2014), have previously consid¬ 
ered radiation from a similar physical setup.® During the time 
of writing also Gao et al. (2015) have investigated EM emis¬ 
sion in a similar context. All of the aforementioned authors 
considered a setup that is qualitatively similar to the phe¬ 
nomenology in Phase III of our model; a remnant NS (pre¬ 
viously born in a BNS merger), surrounded by a PWN that is 
confined by a shell of matter. They present different models 
to evolve such a setup and to predict thermal and non-thermal 
radiation emerging from such a system. Based on a simple dy¬ 
namical model, Yu et al. (2013) find thermal emission peak¬ 
ing around ^ — 10® s after merger with luminosities of 

^ 10"^'*—10^® erg s“^, which they term “merger-nova”. Build¬ 
ing on this model, Gao et al. (2015) argue that such a merger- 
nova is consistent with the late-time rebrightening observed 
at optical and X-ray wavelengths in GRB 080503. Metzger & 
Piro (2014) consider a more detailed physical model, in par¬ 
ticular, regarding the PWN physics, but without relativistic 
dynamics. Eurthermore, they implement a detailed formalism 
to compute the ejecta opacities, the degree of ionization, and 
the resulting albedo in a self-consistent way. As a result of a 
more detailed physical description of the PWN and its inter¬ 
action with the ejecta layer, Metzger & Piro (2014) obtain a 
dimmer thermal optical/UV transient peaking at a luminosity 
of ^ 10"*® — 10"^'* erg s“® on a timescale of ~ 10^ — 10® s af¬ 
ter merger and argue that the late-time X-ray excess of GRB 
130603B (Fong et al. 2014) as well as the late-time optical re¬ 
brightening and X-ray emission of GRB 080503 is consistent 
with their model. 

In contrast to previous work, we start the evolution shortly 
after the BNS merger and introduce a first baryonic wind 
phase that we expect to be generic of BNS mergers leading 
to the formation of a long-lived remnant NS. Furthermore, 
we introduce the pulsar ignition and pulsar wind shock phase 
(qualitatively similar to the initial phase of Metzger et al. 
2014). As the pulsar wind shock is propagating at relativis¬ 
tic speeds in our setup, we implement a relativistic scheme to 
describe the propagation across the ejecta and a detailed de¬ 
scription of the energy transfer between the PWN, the shock 
heated, and unshocked ejecta during this phase. In Phase III, 
we implement a self-consistent model for the radiative pro¬ 
cesses occurring in the PWN, including pair creation and an¬ 
nihilation, Compton scattering, Thomson scattering, and syn¬ 
chrotron cooling. Moreover, we employ a relativistic descrip¬ 
tion of the dynamics in terms of a Milne-universe model and 
develop a reconstruction of the observer lightcurve and spec¬ 
tra taking into account the combined effects of relativistic 
beaming, the relativistic Doppler effect, and the time-of-flight 
effect. As a result of baryon pollution being a general feature 
and the level of detail reached here, we expect our model to 
be applicable to a much larger class of SGRBs than previously 
thought. We explore this possibility in Paper II. In particular, 
we point out that previous magnetar models, which are ap- 

* Gao et al. (2013) have also considered a similar scenario, but focus on 
computing afterglow radiation from the interaction of ejected mass with the 
ambient medium. 


plied to large classes of SGRB events (e.g., Rowlinson et al. 
2013; Gompertz et al. 2013, 2014; Lu et al. 2015) neglect the 
effects of surrounding ejecta material on the magnetar emis¬ 
sion and instead assume a direct and instantaneous conver¬ 
sion of spin-down energy into observed X-ray luminosity 
Lx by some unspecified process, Lsd cx: Lx- This typically 
leads to very simple analytical fitting formulae for the X-ray 
lightcurves, which is in sharp contrast to the self-consistent 
dynamical evolution considered here. 

Rayleigh-Taylor instability, kick velocities, and PWN jet — Our 
model is one dimensional and highly idealized. As noted by 
Metzger & Piro (2014), at early times the high pressure of the 
PWN pushing against the ejecta shell can trigger a Rayleigh- 
Taylor instability. This would presumably cause some poros¬ 
ity in the ejecta envelope thanks to which non-thermal radi¬ 
ation from the PWN could escape from the system already 
at much earlier times. Such effects, however, are difficult to 
take into account in a one dimensional model like ours and 
are neglected here. 

Another deviation from spherical symmetry not accounted 
for by our one-dimensional model could arise from the pos¬ 
sible presence of a kick velocity of the newly-formed NS. A 
large kick velocity could qualitatively alter the evolutionary 
scenario considered here. 

Furthermore, based on axisymmetric magnetohydrody¬ 
namic simulations, Bucciantini et al. (2012) show that for 
sufficiently high spin-down energies, the built-up of strong 
toroidal magnetic field in the PWN interior can drive a bipolar 
jet through the ejecta shell and might even disrupt it entirely. 
Such jet breakouts were employed by Bucciantini et al. (2012) 
to model SGRBs with extended emission. However, recent 
three-dimensional simulations indicate that such jet breakouts 
are an intrinsic effect of two-dimensional simulations, and 
that the kink instability in three dimensions destroys the po¬ 
lar jets and dissolves them into the PWN (Porth et al. 2014). 
Therefore, we exclude such catastrophic events for our model. 

Future improvements — Our model already reflects a certain 
degree of detail and sophistication, but many aspects can be 
improved in future work. A more accurate description of ra¬ 
diative transfer through the ejecta shell and the PWN is re¬ 
quired to predict more accurate lightcurves and spectra, which 
would be desirable for a more detailed comparison with ob¬ 
servational data. Furthermore, the present implementation of 
a self-consistent modeling of the radiative processes occurring 
in the PWN requires the assumption of quasi-stationarity (as 
far as those processes are concerned). In the future, it would 
be desirable to develop a time-dependent formalism, which 
would then allow us to include a self-consistent computation 
of a time and frequency-dependent albedo of the ejecta mate¬ 
rial (similar to Metzger & Piro 2014). Moreover, such a time- 
dependent formalism would more accurately describe further 
acceleration of the ejecta shell in Phase III, it would allow us 
to include the associated pAV work done by the nebula, and 
it would provide a more accurate model of the PWN emis¬ 
sion during the transient phase following the collapse of the 
remnant NS. Finally, we note that with different initial data 
and some modifications, our model could also be employed 
to investigate EM afterglows of long GRBs. 
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APPENDIX 


RELATIVISTIC RANKINE-HUGONIOT CONDITIONS 

This appendix is devoted to deriving Equations (58), (59), (60), and (64) from the special-relativistic Rankine-Hugoniot condi¬ 
tions (e.g., Taub 1948; Marti et al. 1994), thereby also rederiving some of the expressions in Taub (1948) using our notation and 
conventions. The Rankine-Hugoniot conditions were obtained in the special-relativistic case by Taub (1948) and relate the fluid 
properties across a shock wave imposing continuity of the mass and energy-momentum flux: 

= 0, (Al) 

= 0. (A2) 


Here, [/] = /r — /l relates the values of a quantity / on one side and the other side of the discontinuity surface with normal vector 
n^. We choose a frame in which the shock is at rest and adopt Cartesian coordinates such that = (0,1,0, 0). Furthermore, 
we assume an ideal fluid with rest-mass density p, pressure p, specific internal energy e, four-velocity = ( 7 , ju, 0 , 0 ), where 
7 = (1 — and energy-momentum tensor = p(c^ + e + p/p)u^u'' -f where — diag(— 1 , 1 , 1 , 1 ) is the 

Minkowski metric. Finally, we assume an ideal gas equation of state, p = (F — l)pe, where F denotes the adiabatic index. Under 
these assumptions, the Rankine-Hugoniot conditions (Equations (Al) and (A2)) are written as (cf also Taub 1948) 
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We note that for a mixture of radiation and ideal gas with F = 4/3 (i.e., the ejecta matter considered in our model), the Rankine- 
Hugoniot conditions take exactly the same form, with pl and pr being replaced by the respective sums of the radiation and fluid 
pressures. 

Squaring Equation (A4) and subtracting Equation (A5) multiplied by [... -f ...] yields (where [... -f ...] denotes the square 
bracket on the right hand side of Equation (A5) with a + sign separating the two terms; cf. also Equation (7.8) in Taub 1948): 
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(A 6 ) 


This equation can be read as a quadratic equation for pL in terms of pr, pl, and pr. In the following, we assume a strong shock, 
i.e.. 


Pl > Pr, 


(A7) 


and a non-relativistic fluid ahead of the shock, i.e., 

Pr/(Prc^) < 1- (A8) 

These assumptions are very well satisfied for the pulsar wind shock in Phase II (see Section 4.2.3), given the non-relativistic 
unshocked ejecta ahead of the shock front and the very high nebula pressure p„ (Equation 53) that equals the pressure pl of 
the shocked ejecta according to the pressure balance condition (61). With these assumptions, we obtain Equation (60) from 
Equation (A 6 ): 
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With the help of Equation (A3), Equation (A5) can be written as 
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from which using Equations (A7) and (A8) we obtain: 
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In our setup, ur < 0, i.e., the fluid ahead of the shock is moving toward the shock in the rest frame of the shock front. As seen 
from the frame comoving with the fluid ahead of the shock, the shock velocity Ush.r is then given by Ush,R = |wr |c (Equation (59)). 
In order to obtain the shock speed in the lab frame, we apply a Lorentz transformation and arrive at Equation (58). Accordingly, 
we have ul > 0 and the velocity of the fluid behind the shock as seen from the frame comoving with the shock front is given by 
wl = ulc (this quantity is needed in Equation (68)). 

Eor the jump in specific internal energy across the shock, we have 
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Using Equations (A7), (A8), and (A9) this results in 
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which is the desired expression in Equation (64). 


TRANSFORMATION BETWEEN COMOVING AND LAB FRAME 

This appendix defines the lab and comoving frames used in Phase III of our evolution model to describe the expansion of the 
ejecta shell (cf. Section 4.3.2) and it provides transformations for some quantities required by our model. The Milne universe 
metric has recently been employed to describe an expanding homogeneous GRB fireball (Li 2007, 2013). We follow this approach 
and apply it to the expansion of a thin ejecta shell. Some of the expressions that we derive here are based on expressions discussed 
by Li (2007) and Li (2013), some of which we shall rederive here for completeness. 

Let = {ct, r, 9, (p) denote the rest frame of the NS (or the “lab frame”) with spherical coordinates and Minkowski metric 
9p,v = diag(—1,1) sin^ 6). Now consider the following coordinate transformation X''^ = {cq, 9, p), defined by 

f = 77 C 0 sh^, r = cqsinh^, (Bl) 

which transforms the Minkowski metric into = diag(—1, a^{q), a?{q) sinh^ a?{q) sinh^ ^sin^ 9), with scale factor a{q) = 
cq, which is known as the metric of the Milne universe (cf., e.g.. Equation (16.15) of Rindler 2006). Eor a test particle with 
constant radial velocity v, we have r = cpt, where f3 = vjc. Using Equation (Bl) this yields 

^ = arctanh/?, q = (B2) 

where 7 = 1/ is the Lorentz factor. Hence, q is the proper time of the particle and ^ its rapidity. At any time t, the ejecta 

shell can be assigned a Milne universe with velocity u = at r = i?ej by appropriately rescaling the lab frame time coordinate. 
Eor the ejecta shell thickness is typically sufficiently small, i.e., Agj ^ i?ej, 31 any given time we can describe the ejecta as a thin 
shell in this Milne universe, which is what we call the comoving frame at time t. 

Eirst, we determine the three-acceleration a in the lab frame in terms of the corresponding quantity a in the comoving 
frame (cf. Equation (95)). The four-acceleration of a particle in the comoving frame is given by A'^ = d'^X'^/dq'^ = 

{q'^av'jc^ + od)), where w'® = {v'^,v'^,v'^) is the three-velocity, a* = dv''^/dq, v' = and 

a = dv'/dq. Eor the fluid in the ejecta shell, v' = u'* = a® = = 0. Hence, 

A"^ = (0,a«,0,0), 
v' = cqv'^, 
a = cqa^. 

Analogously, in the lab frame the four-acceleration is given by avv^ j(? -f a*)), where u* 

the three-velocity, a® = du®/df, v = gijV^vi, and a = dv/dt. Using = 0, we obtain 

A^ = {q^av/c,q^a,0,0). (B 6 ) 

Moreover, using Equations (B3), (B5), (Bl), and (B2), we have A^ = {dX'^/dX'^)A'^ = {dr/dP)a^ = 7 a and thus we obtain 
the desired expression (cf. Equation (95)) 
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Second, we determine the transformation of area, volume, and internal energy. Henceforth, “com” refers to the comoving 
frame. The surface area of spheres of radius r and corresponding rapidity ^ centered around r = ^ = 0 are invariant under the 
transformation between lab and comoving frames: 


= (?rf' sinh^ ^ 


sin 0 d 0 d 0 = sinh^ ^ = 47 rr^ = Sr- 


(B 8 ) 


The corresponding spherical volumes transform as follows: 

14om(0 = J sin ^ d0d(/) J sinh^ ^ d^ 

/•2« 

= TTc^p^ / (cosh z — 1) dz = 7rc^77^(sinh 2^ — 2^) 

Jo 

= = CV{r), 

where 

3 sinh 2^ — 2^ 3 ^‘^j3 — arctanh /3 

4 sinh^^ 2 7^/33 

Note that our definition of C differs from defined in Li (2013) by a factor 7 “^ = (cosh^)“^. For the volume of the ejecta 
shell we thus find: 

V;j,com = Kom(Cej) " 'l4om(Cn) = C'i^(^ej) " {R^) = CKj- (Bll) 

In order to compute the transformation of total energy of a fluid, let T^i, = (e + j))u^vy jc? + denote the part of the 
energy momentum tensor excluding kinetic and rest-mass energy, with e being the internal energy density, p being the pressure, 
being the four-velocity, and being the metric of flat spacetime. The total energy density as seen by a normal (Eulerian) 
observer with four-velocity perpendicular to the hypersurfaces of constant time t in the Minkowski spacetime is given by 
Ctot = T^yU^-nS = 7 ^e + ( 7 ^ — l)p. For a mixture of ideal gas and radiation, we can rewrite Ctot as 

etot = ^( 47 ^ - l)er + [ 7 ^ + ( 7 ^ - l)(r - l)]ef = ^( 47 ^ - l)e, (B12) 

where and Cf refer to the internal energy density of the radiation held and the gas, respectively. In the second equality in Equa¬ 
tion (B12) we have assumed an ideal gas with adiabatic index E = 4/3 (as for the ejecta material we consider). Consequently, 
the total energy of the fluid (apart from kinetic end rest-mass energy) in a spherical volume as measured in the comoving frame 
is given by (cf. Equation (B9b)) 

-Ecom = etot,com 4 / 0 ™ = 7rc^7^(sinh 2^ - 2^)e, (B13) 


(B9a) 

(B9b) 

(B9c) 

(BIO) 


whereas in the lab frame it is given by = 47 r J etotr^ dr, which after some manipulations can be rewritten as 

E = sinh^ ^ e. (B14) 

Hence (cf. Equations (BIO) and (B9)), 

-^com ^ ^om 1 c\ 

Thanks to Equation (Bll), this transformation also holds for shell volumes. 

Finally, we discuss how luminosities transform between the lab and the comoving frame. Let L denote the luminosity of thermal 
radiation originating from the surface of the ejecta shell. Then using Equation (B 8 ), Lcom = ‘5'Cej^^eff com ~ where 

Teff.com — (16/3)T^onj/(ATcom + 1) IS the effective temperature at optical depth Arcom in the comoving frame. The energy loss 
AScom associated with this luminosity as measured in the comoving frame during a time Ag is given by AE^om = LcomAp = 
Lcom'l~^ At, where we have used Equation (B2) and 7 denotes the Lorentz factor of the ejecta shell. The corresponding energy 
loss as measured in the lab frame is then obtained as (cf. Equation (B15)) LAt = AE = (/“^Ai^coni = Hence, 

we arrive at 

L = r^ 7 ”^icom, (B16) 


which is the desired result to motivate Equations (100) and (101). 


SYNCHROTRON EMISSION 

This appendix derives the relevant expressions needed to include effects of synchrotron emission into our model of the PWN in 
Phase III (Equations ( 86 ), (87), and ( 88 )). The spontaneously emitted spectral power at frequency of a single particle of charge 
q, mass m, and Lorentz factor 7 , accelerated by a magnetic held of strength B, is given by (Rybicki & Lightman 2004; Crusius 
& Schlickeiser 1986) 

p{iy,9,j) = „ sin^— [ iLsjzjdz, (Cl) 

mC^ Jyjy, ^ 
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where c is the speed of light, 0 is the angle between the magnetic field and the direction of particle motion, K& denotes the 

modified Bessel function of order 5/3, and = ( 3 / 47 r)(q_B 7 ^/mc) sin0. Assuming an isotropic distribution of 

particle velocities, we can average over all possible angles for a given speed 7 to obtain the average energy loss rate for a particle 
(cf. Crusius & Schlickeiser 1986): 


= ^ / d6» sin6»p(i^,6»,7) =- 

Jo Jo 


/•47r 




where 


n{z) = W - Wi^iiz)W_i^iiz)]. 


(C2) 
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Here, Wx,fi{z) = z^^^U{Q.h -f /i — A, 1 -f 2/i, z) is the Whittaker function and U denotes the confluent hypergeometric 

function of second kind (Abramowitz & Stegun 1972, Section 13.1). Substituting 1 / with the dimensionless energy x = hv jmc?, 
with h being the Planck constant, one has the relation p{x,^) = {mc^/h)p{iy,j)-, accordingly, we define x^ = hvcirai?. 
Therefore we can write the total cooling rate 7 of a single particle as 


1 f , X , f , 

7(7) =-T = - — -^ / TZ{x/xc)dx. 

mc‘^ J hmc^ J 


Furthermore, the number density of photons emitted per unit time and per unit normalized energy x is given by 

VSq^B 1 


n{x) =[ N{'j)p{x, 7 ) d 7 = , f - f N{-f)TZ{x/xc)d'Y 
mc^x J nmc^ x J 


(C4) 
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where N ( 7 ) is the particle density per normalized energy 7 . Finally, assuming a uniform opacity coefficient the optical depth Ar 
to synchrotron self-absorption for a spherical volume of radius R is given by (Ghisellini et al. 1988; Ghisellini & Svensson 1991) 


Ar(i^) = 


R f A( 7 ) d 




7 a dt 


[yap{iy, 7 )] d 7 , 


where a = — Using the information from above. Equation (C 6 ) can be rearranged to give 


At(x) = 


^q^h^BR 1 
x'^ 


N{l) 


d_ 

dj 


TZ{x/x^) + f{'y)Tl{x/x^) 
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where f{-y) = ( 27 ^ - l)/[ 7 ( 7 ^ - 1 )], with lim^_>i /( 7 ) = 2 . 


ENERGY CONSERVATION IN THE NEBULA 

This appendix demonstrates that Equations (78) and (79) conserve energy. Integrating Equation (78) the total energy balance 
per unit volume per unit time is given by 

From Equation (83) we have 

/ 'Tmax rimay. P Plmax 

7c.syn(7)^(7)d7= J N {-f)TZ{x / Xc) - J 7cA^(7)d7 


^'7nitix 


= / nsynX dx — 


icN{i) d7, 


and hence 
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ri^ynX dx = 


[< 3 ( 7 ) + -P(7)](7- l)d7' 




7cW(7) d7. 


(D2) 

(D3) 

(D4) 


Plugging this into Equation (Dl) and noting that the second term on the right hand side of Equation (D4) cancels the third and 
sixth term on the right hand side of Equation (Dl) (see LZ87, Appendix A), we can proceed as in Appendix A of LZ87 to 
conclude that 


TT-escS: dx= rioX dx + 


Q(7)(7- l)d7. 


(D5) 


Therefore, the total injected power equals the total power output, which shows that energy is conserved in steady state. 














